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1BST1ACT 

Sine*  widespread  use  of  the  finite  element  aethed  began 
^  in  the  early  1960*s,  each  effort  has  been  devoted  to  the 

development  of  the  aethod  itself,  while  only  recently  has 
there  been  any  research  directed  at  aininizing  the  discreti¬ 
zation  error  by  a  proper  selection  of  the  eleaent  grid.  This 
paper  ezaaines  soae  recently  proposed  grid  optiaization 
techniques  and  applies  then  to  soae  one-  and  two-diaensional 
linear  self-adjoint  boundary  value  probleas.  Guidelines 
requiring  siniaal  software  aodifisation  are  recoaaended  to 
assist  the  analyst  in  obtaining  iaproved  finite  eleaent 
solutions. 
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I.  IMT8QP0CTI3M 

The  critical  concern  in  an?  finite  eleaent  analysis  is 
the  adequacy  of  the  selected  aesh  to  provide  reliable  solu¬ 
tion  results  within  soae  reasonable  cost  range.  The  goal  of 
finite  eleaent  grid  optiaization  then  becoaes  one  of 
obtaining  aaxiaua  solution  accuracy  for  a  prescribed  anal¬ 
ysis  cost,  while  this  objective  is  generally  not  realized  in 
today's  widespread  use  of  finite  eleaent  analysis,  the  effi¬ 
cient  coaputation  of  solutions  with  optimal  accuracy  will 
becoae  paraaount  to  the  engineer  as  finite  eleaent  nethcds 
are  applied  to  increasingly  difficult  dynamic,  nonlinear, 
and  evolution  probleas. 

1.  HISTOBIC1L  B1CKGBOOHD 

In  the  early  1960's,  with  the  help  of  the  high  speed 
digital  coaputer,  finite  eleaent  aethods  revolutionized 
problea  solving  in  engineering.  Since  that  tine  the  major 
research  efforts  have  concentrated  on  expanding  the  theoret¬ 
ical  basis  of  the  aethod  and  extending  its  application  in  a 
variety  of  fields.  Only  recently  has  there  been  significant 
attention  directed  at  niniaizing  finite  eleaent  solution 
errors  by  properly  defining  the  eleaent  grid.  Early  attempts 
at  distributing  the  nodes  and  choosing  the  elements  to 
ensure  sene  degree  of  confidence  in  the  solution  accuracy 
were  predoainantly  dependent  upon  the  analyst's  engineering 
judgeaent  and  experience,  since  there  were  no  established 
procedures  foe  accoaplishi ng  this  objective.  Even  these 
attempts  towards  grid  optiaization  have  become  largely 
ignored  with  the  advent  of  automatic  mesh  generators,  which 
have  drastically  reduced  preprocessing  costs  while 


acco»plishing  little  in  improving  solution  accuracy.  These 
prograss  automatically  construct  the  element  grid,  usually 
in  a  uniform  manner,  after  the  user  merely  defines  the 
problem  and  specifies  the  number  of  elements  to  be  used.  To 
establish  convergence  and  achieve  the  desired  solution  accu¬ 
racy,  the  user  simply  repeats  the  analysis  using  a  finer 
mesh  of  uniformly  distributed  elements  while  monitoring  such 
convergence  indicators  as  successive  solution  values  at 
common  nodal  points  or  the  asymptotic  increase  in  the  energy 
content  of  the  mesh.  The  often  disastrous  flaw  in  such  a 
practice  is  that  for  nearly  degenerate  problems  which 
exhibit  very  large  gradients,  the  asymptotic  range  is  only 
entered  for  an  extremely  large  number  of  degrees-of-freedom, 
often  beyond  computer  limitations  [Bef.  1].  In  this  case, 
uniform  mesh  refinement  may  produce  apparent  convergence, 
when  in  fact  the  solution  accuracy  is  poor.  Therefore,  the 
need  for  a  finite  element  grid  optimization  procedure  is 
clearly  evident. 

The  first  formal  attempts  at  finite  element  grid  optimi¬ 
zation  did  not  begin  until  the  early  1970 *s.  These  early 
approaches  involved  the  inclusion  of  the  nodal  coordinates 
as  dependent  variables  in  the  minimization  of  the  potential 
energy  functional  [Bef.  2].  Unfortunately,  the  resulting 
system  of  equations  is  highly  nonlinear  and  the  computa¬ 
tional  effort  involved  in  its  solution  is  so  great  that 
similar  accuracy  can  be  obtained  at  a  fraction  of  the 
expense,  simply  by  employing  a  very  fine  mesh.  Clearly,  this 
method  does  not  approach  the  finite  element  grid  optimiza¬ 
tion  goal  of  achieving  a  specified  solution  accuracy  for  a 
minimum  analysis  cost.  For  this  reason,  virtually  all  of  the 
grid  optimization  techniques  since  developed  are  based  on  a 
"near -optimum"  strategy  whereby  nearly-optimal  solution 
results  are  obtained  without  the  computational  inefficiency 
of  a  numerical  optimization  analysis.  The  growing  emphasis 
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has  been  on  adaptivity,  a  procedure  for  efficient  construc¬ 
tion  of  near-optimum  grids  by  the  iterative  application  of 
sose  criterion,  based  on  data  already  computed  from  the 
solution  for  a  previous  grid.  This  procedure  is  far  more 
efficient  than  the  conventional  approach  of  repeating  the 
analysis  using  successively  finer  uniform  grids. 
Experimental  self-adaptive  finite  element  codes  have 
recently  been  developed,  starting  from  the  user’s  initial 
idealization,  these  programs  automatically  generate  a  near- 
optimum  grid  and  solve  the  resulting  equations. 


B.  IiVESTIGATITB  APPHOACH 

In  undertaking  any  numerical  optimization  task,  the 
analyst  must  first  define  the  objective  along  with  any 
constraints  to  be  imposed  upon  the  objective  variables;  and 
finally  a  numerical  algorithm  must  be  prescribed  to  perform 
the  optimization,  preferably  one  which  will  do  so  effi¬ 
ciently  for  the  particular  problem.  Since  the  term  "optimum" 
most  often  refers  to  a  solution  obtained  by  mathematical 
programming,  which  is  very  inefficient  for  grid  optimiza¬ 
tion,  a  near-optimum  grid  obtained  by  application  of  an 
adaptive  procedure,  henceforth  will  be  termed  an  optimum 
grid.  However,  before  such  a  grid  can  be  determined  to 
satisfy  the  stated  objective  of  obtaining  maximum  solution 
accuracy  for  a  prescribed  cost,  the  terms  "accuracy"  and 
"cost"  must  be  defined;  but,  more  importantly,  the  optimiza¬ 
tion  goals  must  be  specified.  This  is  critical  because  grid 
optimization  can  be  implemented  in  various  forms  depending 
upon  the  optimization  goals,  which  will,  in  general,  be 
determined  by  the  original  purpose  for  performing  the  finite 
element  analysis  [Bef.  3].  For  example,  if  the  purpose  of 
the  analysis  is  to  evaluate  a  local  quantity,  such  as  the 
maximum  value  of  the  dependent  variable  or  one  of  its 


derivatives,  then  the  nodal  distribution  should  be  densest 
in  the  region  where  that  aaxiaua  is  determined.  If,  on  the 
other  hand,  the  interest  is  on  an  integral  quantity  of  the 
dependent  variable  ever  a  region  of  the  doaain,  then  the 
nodes  should  be  assigned  to  achieve  a  unifora  distribution 
of  the  error  over  that  region.  For  the  purpose  of  this 
investigation,  attention  will  be  focused  on  the  three  finite 
eleaent  resultants  with  the  aost  significance  in  solid 
aechanics  and  other  fields  in  which  finite  eleaent  analysis 
is  eaployed:  the  aaxiaua  value  of  the  dependent  variable,  or 
solution;  the  aaxiaua  value  of  the  gradient  of  the  solution; 
and  the  integral  of  the  solution  over  the  doaain. 

In  order  to  define  the  solution  accuracy,  it  will  be 
necessary  to  coapare  the  error  in  the  solution  resultant 
obtained  using  an  optiaal  grid  to  the  error  obtained  using 
soae  baseline  grid  with  the  sane  number  of  degrees- of- 
freedom.  For  convenience,  the  reference  grid  chosen  will  be 
a  uniform  grid,  or  one  with  all  elements  of  the  same  erder 
and  approximately  the  same  size,  with  the  understanding  that 
no  knowledgeable  analyst  would  attempt  to  use  such  a  grid  in 
the  solution  of  a  problem  with  large  gradients. 

The  definition  of  analysis  cost  will  be  greatly  simpli¬ 
fied  by  making  the  assumption  that  it  is  directly  propor¬ 
tional  to  the  number  of  degrees-of- freedom  used  to  obtain 
the  solution.  In  reality  the  cost  depends  on  many  factors, 
soae  of  which  are  very  difficult  to  quantify. 
Understandably,  the  number  of  degrees-of-freedom  is  not  the 
sole  measure  of  computational  costs,  but  it  appears  to  be  an 
adequate  aeasure  of  preprocessing  and  postprocessing  costs 
which  often  account  for  the  major  portion  of  the  total 
analysis. 

This  investigation  will  concentrate  on  the  use  of  finite 
eleaent  grid  optimization  methods  for  solving  problems  of 
structural  mechanics.  While  this  area  has  dominated  the 
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literature  on  the  subject,  the  techniques  presented  herein 
extend  equally  as  veil  to  any  field  for  vhich  variational 
principles  apply. 

There  are  two  hey  questions  vhich  aust  be  ansvered  prior 
to  the  adaptive  application  of  finite  element  grid 
optiaization: 

(1)  What  criterion,  based  on  the  results  of  an  initial 
finite  eleaent  analysis,  should  be  used  to  indicate 
regions  vhere  the  original  grid  is  inadequate  ? 

(2)  Hhat  method  of  grid  refinement  should  be  employed  ? 

Considerable  attention  vill  be  devoted  to  these  questions  in 
the  next  tvo  chapters. 

C.  OBJECTIVES 

The  objectives  of  this  investigation  are: 

(1)  To  examine  sone  recently  developed  grid  optiaization 
techniques  vhich  have  reached  the  state  of  practical 
applicat  io  n. 

(2)  To  implement  some  of  these  techniques  in  the  solu- 
tion  of  some  one-  and  tvo-diaensional  linear  self- 
adjoint  boundary-value  problems. 

(3)  To  drav  a  comparison  among  these  applied  techniques 
in  terms  of  solution  accuracy,  analysis  cost,  and 
ease  of  implementation. 

(4)  To  recommend  sone  guidelines  to  assist  the  analyst 
in  obtaining  optimal  finite  element  solutions 
employing  currently  available  or  easily  amendable 
so ft vara. 
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The  primary  theoretical  concern  in  the  application  of 
adaptive  grid  optimization  is  the  selection  of  the  refine¬ 
ment  criterion.  In  ether  words,  one  must  decide  upon  which 
solution  parameter,  obtained  from  an  initial  idealization, 
nay  most  appropriately  be  used  to  give  some  indication  as  to 
where  the  initial  grid  is  inadequate  and  thus  needs  refine¬ 
ment.  There  are  several  competing  proposals  concerning  the 
most  appropriate  choice  of  a  refinement  criterion.  In 
reality,  the  decision  must  be  based  upon  such  factors  as  the 
type  of  problem  being  solved,  which  criterion  is  most  prac¬ 
tically  implemented,  and  whether  accuracy  is  desired 
locally,  globally,  pointwise,  or  with  respect  to  an  integral 
quantity.  The  following  are  some  of  the  more  practical 
refinement  criteria  used  in  grid  optimization  at  present. 

A.  SOLUTION  PABAHBTBB  fABIATIO* 

The  most  direct,  computationally  inexpensive,  and 
earliest  proposed  indication  of  where  an  element  grid 
requires  refinement  is  a  measure  of  the  variation  of  some 
solution  parameter  over  the  domain.  It  is  only  logical  that 
a  piecewise  polynomial  approximation  would  experience  the 
most  difficulty  in  modeling  the  desired  response  in  those 
regions  where  the  solution  or  its  resultants  were  either  not 
smooth  or  were  characterized  by  large  gradients.  Therefore, 
the  basis  of  this  criterion  is  to  refine  the  mesh  in  those 
areas  where  a  solution  parameter  varies  rapidly,  with  the 
implication  that  the  optimum  mesh  is  one  for  which  the  solu¬ 
tion  parameter  variation  over  each  element  is  uniform 
throughout  the  domain.  The  first  consideration  in  the 
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application  of  such  a  'criterion  is  to  find  a  scheme  for 
distributing  the  nodes  to  achieve  such  a  condition.  For 
one- dimensional  problems  the  task  is  trivial,  but  one  way  to 
ensure  equal  variation  over  each  eleaent  in  higher  dimen¬ 
sions  is  to  position  the  nodes  along  equidistant  contours  of 
the  chosen  paraseter.  This  is  precisely  the  procedure 

prescribed  for  a  practical  optimization  technique  known  as 
contouring.  The  other  consideration  is  the  determination  of 
which  solution  parameter  is  to  be  used.  In  fact,  several 
solution  parameters  have  been  found  to  work  quite  well 
[Hef.  4],  but  the  most  commonly  used  and  the  one  that  is 
consistent  with  the  potential  energy  minimization  formula¬ 
tion  is  the  strain  energy  density  [  Hef .  2].  Because  its 
employment  requires  only  minor  software  changes  and  it  has 
been  found  to  produce  excellent  results,  this  refinement 
criterion  was  used  extensively  throughout  the  course  of  this 
research. 

B.  GRID  ITEH1TIOI 

Another,  rather  basic  but  less  direct,  method  of 
locating  regions  of  the  mesh  to  be  refined  is  known  as  grid 
iteration,  which  can  be  implemented  in  one  of  two  ways.  An 
initial  coarse  grid  analysis  may  be  repeated  using  either  a 
finer  or  a  higher  order  mesh,  and  a  comparison  of  the  resul¬ 
tants  of  interest  between  the  two  solutions  will  identify 
those  areas  of  the  domain  where  refinement  is  most  effec¬ 
tive.  Another  approach  is  based  on  the  assumption  that  the 
greatest  benefit  is  to  be  gained  by  refinement  in  those 
regions  where  a  small  perturbation,  like  the  introduction  of 
a  single  additional  degree- of-freeiom ,  produces  the  greatest 
change  in  the  solution  or  one  of  its  resultant  parameters. 
Since  additional  degrees-of-f reedom  would  be  expected  to 
produce  the  greatest  change  in  those  regions  where  the 


desired  response  varied  most  rapidly,  refinements  based  or. 
this  method  provide  results  very  similar  to  those  obtained 
using  the  solution  parameter  variation  criterion  already 
discussed.  The  grid  iteration  method  nay  at  first  appear  to 
be  more  computationally  expensive,  but  it  was  developed  to 
be  most  efficiently  implemented  employing  a  special  family 
of  elements.  These  elements,  called  "hierarchical",  possess 
some  very  desirable  properties  for  this  application  and  will 
be  discussed  in  the  next  chapter. 

C.  BLBBB1T  BESIDOALS 

The  major  drawback  with  refinement  criteria  based  on 
solution  parameter  variations  is  that  their  applicability 
appears  limited  to  elastostatic  problems.  For  this  reason, 
several  investigators  have  recently  developed  grid  refine-* 
neat  criteria  based  on  element  residaals,  which  appear  prem¬ 
ising  for  application  to  all  types  of  finite  element 
problems,  including  nonlinear  analysis.  The  reason  for  this 
is  that  the  rasidual  has  essentially  the  same  meaning 
regardless  of  the  problem  type  [Bef.  5].  For  example, 
consider  the  governing  differential  equation, 

D  (  u  ]-f»0 

defined  cn  some  domain,  where  3  [  ]  is  a  linear  or  nonlinear 
differential  operator,  and  the  dependent  variable  u  and  the 
non- homogeneous  term  f  are  both  functions  of  the  independent 
variables.  let  the  finite  element  approximation  to  the 
solution  of  the  differential  equation  be  ^  ^  u.  Applying  the 
differential  operator  to  the  approximation  gives  rise  to  the 
residual,  which  is  defined  as 

b  *  o  {  a  ]  -  f 

and  is  not  identically  zero  unless  the  finite  element  solu¬ 
tion  is  exact. 
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The  key  to  using  the  residual  as  a  criterion  for  deter* 
■ining  regions  of  the  doaain  where  refineaent  is  necessary 
is  the  local  residual  on  the  eleaent  level,  which  indicates 
the  contribution  of  the  eleaent  to  the  total  error  of  the 
finite  eleaent  approxiaatinn.  Since  the  residual  is  a  point- 
wise  quantity,  the  useful  aeasure  of  the  eleaent  error 
contribution  is  a  nora  of  the  eleaent  residual,  or  the  inte¬ 
gral  over  the  eleaent  of  the  product  of  the  residual  and 
soae  weight  function.  The  integration  is  aost  readily 
perforaed  using  nuaerical  quadrature.  The  grid  optiaization 
strategy  then  becoaes  one  of  refining  the  aesh  so  as  to 
equi -distribute  the  eleaent  residual  noras,  by  forcing  then 
to  becoae  saaller  and  aora  unifora  over  the  doaain. 

There  are  soae  drawbacks  to  the  eleaent  residual  refine- 
aent  criterion  which  have  not  yet  been  fully  resolved,  such 
as  appropriate  selections  of  the  residual  nora  and  the 
weight  function,  and  in  the  conputation  of  the  residual 
itself.  While  the  evaluation  of  the  residual  presents  no 
particular  difficulties  in  the  interior  of  the  eleaent,  it 
is  rarely  deterainable  at  the  edges.  The  reason  for  this  is 
that  in  foraulating  the  finite  eleaent  approxiaation  the 
eleaent  shape  functions  are  defined  so  as  to  provide  only 
that  degree  of  continuity  required  to  adequately  aodel  the 
physical  problea;  the  aost  frequent  consequence  being  that 
D  [ u  ]  is  undefined  along  the  intereleaent  boundaries. 
Unfortunately,  this  singularity  cannot  be  ignored  and  a  more 
coaplicated  analysis  aust  be  applied  in  order  to  bound  the 
residual  contributions  at  these  boundaries  [Bef.  6]. 

0.  1- POSTS *10 B I  B1B0B  BSTZHiTBS 

A  sophisticated  extension  o  f  the  eleaent  residual 
criterion  is  one  based  on  coaputable  error  estinates  from  an 
initial  finite  eleaent  analysis.  This  utilizes  the  energy 


17 


\ 


non  of  the  residual,  in  which  case  the  weight  function  is 
the  residual  itself.  Sesearch  in  reliable  error  estimates 
was  pioneered  by  Babuska  [Ref.  7,  8]  for  linear  quadrilat¬ 
eral  elements  and  aore  recently  by  Zienkiewicz  [Ref.  9]  for 
higher  order  eleaents.  The  aajor  difference  froa  the 
residual  criterion  is  that  instead  of  equi-distributing  the 
eleaent  residual  nor as  over  the  doaain,  they  are  noraalized 
to  coapute  error  indicators  for  the  eleaents,  which  are  in 
turn  used  to  coapute  reliable  poiatwise  error  estiaates  for 
the  solution  as  well  as  the  energy  error  over  the  doaain. 
These  quantities  are  of  priaary  iaportance  because  they 
provide  not  only  an  indication  of  where  additional  refine- 
aent  is  aost  effective,  but  also  a  aeasure  of  the  quality  of 
the  aesh  to  deteraine  whether  additional  refineaent  is 
necessary  [Ref.  9].  The  optiaization  strategy  is  to  obtain 
a  nearly  unifora  distribution  of  the  error  indicators 
throughout  the  doaain,  which  corresponds  to  ainiaizing  the 
error  in  the  energy  nora.  The  refineaent  procedure  aay  prog¬ 
ress  until  all  the  lccal  errors  are  within  soae  prespecified 
tolerance.  Hhile  the  practical  utility  of  such  a  refinement 
criterion  is  obvious,  the  aatheaatical  development  and  the 
algorithas  involved  are  rather  complicated.  However,  the 
process  is  not  computationally  expansive,  and  there  now 
exists  a  prototype  self-adaptive  finite  element  code,  PEARS, 
which  iapleaents  this  refinement  criterion  [Ref.  10]. 


•Vv :  -■> 
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One*  it  has  bssn  determined  where  the  initial  aleaent 
grid  is  inadequate  and  needs  refineaent,  the  next  considera¬ 
tion  is  how  the  idealization  in  these  areas  should  be 
improved.  The  choice  of  the  refineaent  method  to  be  employed 
may  veil  be  a  more  important  decision  than  the  selection  of 
one  of  the  refineaent  criteria  previously  discussed,  since 
at  least  one  investigator  has  observed  that  for  a  particular 
method  of  grid  refineaent,  the  various  refineaent  criteria 
produce  essentially  the  same  solution  results  [Bef.  11]. 

Grid  refineaent  is  the  process  of  introducing  additional 
degrees-of- freedom  into  selected  regions  of  the  finite 
element  grid,  and  may  be  performed  by  one  of  three  methods: 

(1)  The  polynomial  degree  of  the  elements  remains  fixed, 
usually  at  a  lov  value,  while  the  size  of  the 
elements  is  reduced.  This  has  become  known  as  the 
h-version  since  element  size  is  comnonlv  denoted  by 
the  letter  h. 

(2)  The  size  of  the  elements,  usually  few  in  number, 
remains  fixed  while  the  polynomial  degree  of  the 
elements  is  increased.  This  has  become  known  as  the 
p-version  since  polynomial  order  is  commonly  denoted 
by  the  letter  p. 

(3)  The  size  of  the  elements  may  be  reduced  concurrently 
with  an  increase  in  their  polynomial  order.  This  is 
known  as  the  ccabined  h-  and  p-version  of  the  finite 
element  method. 
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I.  CCIVBBG1BCE  OF  6BID  HBFIBEHBIT 


It  is  vail  known  that  the  finite  element  method 
converges  with  an  increasing  number  of  degrees-of- freedom; 
in  fact*  this  is  the  justification  for  its  development, 
therefore*  the  appropriate  aeasure  of  the  effectiveness  of  a 
particular  grid  refineeent  aethod  should  be  the  associated 
rate  of  convergence*  which  generally  will  be  affected  by  the 
seoothness  of  the  approxiaated  function  over  the  subdoaain 
of  interest.  It  has  been  desonstrated  that  when  the  refine- 
aent  is  perfocaed  uniformly  over  the  doaain*  the  associated 
rate  of  convergence  is  asyaptotic*  provided  the  nuaber  of 
degrees-of-freedoa  is  sufficiently  large  [Bef.  1].  The 
asymptotic  rate  of  convergence  is  often  aeasured  as  the 
slope  of  the  error  versus  cost  curve  in  the  linear*  or 
asyaptotic*  range  when  plotted  on  a  log-log  scale.  In 
perforaing  such  an  error  analysis  for  the  displacement 
formulation  of  the  finite  element  method*  the  error  is 
usually  the  relative  strain  energy  error*  approxiaated  by 
the  energy  nora*  and  the  cost  is  assumed  to  be  some  simple 
function  of  the  average  eleaent  size  or  the  nuaber  of 
degrees-of-freedoa  [Bef.  12:  p.  726].  Only  in  the  past 
several  years  has  there  been  any  significant  research 
comparing  the  relative  aerits  of  the  different  methods  of 
grid  refineaent.  Since  the  solutions  of  elliptic  boundary 
value  problems  are  usually  very  smooth  over  convex  domains 
except  in  the  vicinity  of  corners*  aost  of  this  research  has 
focused  on  solutions  exhibiting  singularities*  which 
severely  hinder  the  rate  of  convergence*  as  in  problems  of 
fracture  mechanics  and  in  doaains  with  re-entrant  corners 
[Bef.  1*  13*  14*  15]. 

In  order  for  a  finite  eleaent  analysis  to  be  both  effi¬ 
cient  and  reliable*  the  asyaptotic  convergence  range  should 
be  entered  for  as  few  degrees-of -freedom  as  reasonably 


possible.  In  general  the  p-version  satisfies  this  require- 
sent  better  than  the  h-version.  Sfhile  it  has  been  estab¬ 
lished  that  p-ccnvergence  will  necessarily  occur  whenever 
h-convergence  occurs,  the  converse  is  not  true.  For  example, 
if  the  h-version  using  a  uniform  grid  of  linear  elements  is 
applied  to  a  nearly  degenerate  problem,  the  number  of 
degrees- of- freed  cm  required  for  entry  into  the  asymptotic 
range  may  be  beyond  the  computer's  cound-off  limitations,  in 
which  case  convergence  will  not  occur  unless  the  polynomial 
order  is  increased  [Ref.  1].  Numerical  experiments  on  such 
problems  clearly  indicate  that  the  p-version  requires 
considerably  fewer  degrees- of-freedom  than  the  h- version  to 
achieve  the  same  degree  of  accuracy.  Recent  analyses 
[Ref.  1,  13]  of  the  asymptotic  rata  of  convergence  in  energy 
for  non-smooth  solutions,  using  uniform  refinement  with 
sufficiently  high  numbers  of  degrees-of -freedom,  have  demon¬ 
strated  that  the  p-version  cannot  have  a  slower  rate  of 
convergence  than  the  h-version.  Furthermore,  if  the  singu¬ 
larity  is  confined  to  element  boundaries,  as  is  usually  the 
case,  the  error  for  p-method  is  inversely  proportional  to 
the  number  of  degreem-of- freedom,  whereas  the  error  is 
inversely  proportional  to  the  square  root  of  the  number  of 
degrees-of-freedos  in  the  h-version.  In  other  words,  for 
this  special  class  of  problem,  tne  rate  of  convergence  for 
the  p-version  is  twice  that  for  the  h-version,  which  is  due 
primarily  to  the  ability  of  higher  order  polynomials  to 
"absorb"  singularities  occurring  at  the  element  boundaries. 
This  implies,  at  least  for  this  type  of  problem,  that  in 
order  to  minimize  the  error  for  a  specified  number  of 
degreem-of-freedom,  the  best  strategy  is  not  to  subdivide 
the  domain  uniformly,  but  to  use  instead  a  single  element  of 
increasing  polynomial  order  [Ref.  15]. 
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Sine*  it  is  unlikely  that  one  would  attaapt  to  solve 
such  a  problea  using  uniformly  finac  grids,  a  sore  useful 
coaparison  between  the  convergence  rates  of  the  two  versions 
would  be  based  on  adaptive  refineaent  employing  one  of  the 
solution -based  criteria  discussed  in  the  previous  chapter. 
It  so  happens  that  the  h-version,  when  used  with  optiaally 
refined  meshes,  can  have  a  higher  convergence  rate  than  the 
uniforaly  distributed  p-version,  provided  that  the  element 
order  is  sufficiently  high.  However,  the  p-version  can  also 
be  employed  with  an  optimal  refineaent  criterion.  while 
there  are  yet  no  proven  theorems  concerning  the  convergence 
rates  for  non-uniform  refineaent,  obtaining  the  desired 
solution  accuracy  with  optimal  p-distri  but  ions  appears  to  be 
auch  less  sensitive  to  the  particular  choice  of  the  elements 
to  be  refined  than  with  optimal  h-ref inament  [Ref.  13]. 

It  would  3eea  plausible  that  an  even  better  optimization 
strategy  would  involve  a  proper  combination  of  both  the  h- 
and  p-versions.  It  has  been  demonstrated  for  problems  with 
corner  singularities,  that  a  proper  sequence  of 
h-ref ineaents  combined  concurrently  with  the  proper  sequence 
of  p-distributious  results  in  extremely  high  convergence 
rates,  conjectured  tc  be  exponential  [Ref.  15].  However, 
this  proper  combination  is  difficult  to  determine,  and  adap¬ 
tive  refinement  based  on  the  combined  h-  and  p-versiens 
poses  some  difficult  data  management  problems.  To  avoid  this 
problea  a  more  promising  approach,  proposed  by  Babuska  and 
Szabo  [Ref.  1],  employs  a  graded  mesh  in  which  the  element 
sizes  are  first  reduced  according  to  a  geometric  progression 
towards  the  singularity,  followed  by  determining  the  optimal 
p-distribution  for  those  elements  using  an  adaptive 
criterion.  However,  obtaining  the  optimal  combination  when 
employing  this  scheme  can  be  a  delicate  matter  and,  astcund- 
ingly,  the  highest  accuracy  is  achieved  when  the  polynomial 
order  of  the  elements  actually  increases  with  distance  from 
the  singularity. 
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There  are  soie  additional  advantages  of  the  p-version 
worth  aentioaing.  Because  the  p-version  employs  fewer 
eleaent s,  there  are  lesser  preprocessing  and  postprocessing 
costs  than  for  the  h-version.  Furthermore,  when  bandwidth 
Minimization  and  sparse  Matrix  solution  techniques  are  used, 
the  solution  tiae  for  the  p-version  is  approximately  the 
saae  as  for  the  h-version  for  a  specified  number  of 
degrees-of- freedom,  and  the  p-version  appears  less  suscep¬ 
tible  to  round-off  errors.  Finally,  the  p-version  is  simpler 
to  iaplaaent  adaptively  than  the  h-version  when  hierarchical 
elements  are  employed  [Ref.  13]. 

B.  HIERARCHICAL  FI HITS  BLB HERTS 

The  hierarchical  concept  was  first  introduced  as  a 
simple  method  for  implementing  the  p-version  and  as  a 
convenient  device  fcr  imposing  boundary  continuity  between 
elements  of  different  polynomial  order  [Ref.  9].  Since  then 
a  useful  family  of  elements  based  on  the  hierarchical 
concept  has  baen  developed  and  incorporated  into  COHET-x,  an 
experimental  finite  eleaent  code,  developed  by  Szabo,  which 
self -adaptively  employs  both  the  h-  and  p-versions  of  the 
finite  element  method  [Ref.  IS]. 

For  a  brief  description  of  the  hierarchical  concept 
consider  the  conventional  finite  element  formulation  which 
produces  the  following  system  of  aquations: 

K.  .u(n)*  f<n)  (Eqn.  3.1) 

%(nH  ^ 

where  n  is  the  number  of  degrees-of- freedom,  is  the 

n  x  n  global  stiffness  matrix,  is  the  finite  element 

approximation  of  the  exact  solution,  and  is  the 

n-component  global  load  vector.  Rhen  a  higher  order 


23 


degrees-of-freedcm  are  added  to  the  original  system  using 
conventional  refinement  methods  rhe  system  of  equations 

becomes: 


•^(n-hn) 

% 


where  the  contributions  to  K 


(n-Hn) 


and  f 

% 


(n-Hn) 


(Egr.  3.2) 


from  the  refined 


elements  result  in  a  completely  different  set  of  coeffi¬ 
cients.  If,  on  the  other  hand,  this  refinement  had  been  made 
hierarchically,  the  equations  would  become: 


K,  ,  K,  ,  u(n) 

%(n)  fy(n,m)  ^ 


K  K  UVJ 

%( mfn)  % 


(Eqn.  3.3) 


where  and  j^are  the  stiffness  matrix  and  force  vector 
from  the  original  system  of  equations  for  n  degrees-of- 
freedom  appearing  in  Equation  3.1.  However,  u^m)is  not  the 
nodal  values  of  the  finite  element  solution  for  the  m 
additional  degrees-of-freedom,  but  instead  represents  the 
difference  between  those  values  and  the  pointwise  values 
obtained  from  the  lower  order  polynomial  interpolation  for 
the  unrefined  mesh  of  n  degrees-of-freedom. 

The  primary  advantage  of  hierarchical  elements  is  imme¬ 
diately  observable  from  Equation  3.3.  Because  the  shape 
functions  of  an  element  of  order  p  constitute  a  subset  of 
the  shape  functions  of  an  element  of  order  p+1,  the  local 
stiffness  matrix  and  force  vector  for  each  hierarchical 
element  is  embedded  in  the  stiffness  matrices  and  force 
vectors  of  all  hierarchical  elements  of  higher  order. 

Therefore,  the  global  stiffness  matrix  K/^and  force  vector 

(n)  ^ 

f  of  the  original  system  are  preserved,  thus  saving 


considerable  time  and  effort  expended  on  computing  the  coef¬ 
ficients  for  successive  refinements  [Hef.  14].  Another 
advantage  is  that  the  hierarchical  form  of  the  global  stiff¬ 
ness  matrix  is  more  diagonally  dominant  than  the  one 
resulting  from  a  conventional  refinement,  resulting  in 
improved  conditioning  and  faster  convergence  when  iterative 
solvers  are  employed  [Bef.  9].  Another  benefit  of  hierar¬ 
chical  elements,  which  arises  from  the  "add-on'*  nature  of 
the  nodal  variables  cf  the  higher  order  degrees-of- freedom, 
is  that  the  problem  of  maintaining  boundary  continuity 
between  elements  of  different  polynomial  order  becomes 
trivial.  Instead  of  introducing  global  constraint  equations 
for  the  higher  order  degrees-of-freadom,  the  nodal  variables 
are  simply  set  equal  to  zero  and  condansed  out,  as  if  they 
were  zero-valued  Dirichlet  boundary  conditions  [Bef.  2]. 

There  are  two  major  drawbacks  with  hierarchical  elements 
that  have  hampered  their  widespread  acceptability.  The 
first,  which  has  already  been  mentioned,  is  that  the  nodal 
variables  for  the  higher  order  degrees-of-freedom  represent 
difference  values  rather  than  the  more  easily  identifiable 
values  of  the  dependent  variable  itself.  Secondly,  when 
implementing  the  h-version  of  the  finite  element  method, 
special  integration  rules  must  be  introduced  when  the  subdi¬ 
vided  element  is  in  hierarchical  fora  [Bef.  9].  Of  course, 
the  latter  problem  can  be  evaded  by  using  the  p-version,  for 
which  the  hierarchical  concept  was  developed.  In  spite  of 
the  disadvantages  of  hierarchical  elements,  their  consider¬ 
able  computational  efficiency  and  utility  for  grid  optimiza¬ 
tion  will  certainly  result  in  their  widespread  utilization 
in  future  adaptive  finite  element  nodes. 


i7.  ssxj)  omamnog  issikqsis 


Once  the  analyst  has  identified  where  the  initial  grid 
needs  enrichsent  and  decided  which  refineaent  aethod  to 
employ,  be  aast  then  deter  nine  a  systeaatic  procedure,  or 
algorithm,  to  perfora  the  refineaent  according  to  the 
criterion  selected.  The  ultimate  goal  of  such  a  procedure  is 
to  design  an  element  grid  which  meets  the  optimization 
objective  of  obtaining  maximum  solution  accuracy  for  a  spec¬ 
ified  analysis  cost.  While  the  analyst  may  or  may  not  have 
an  indication  of  the  accuracy  of  the  solution,  he  should 
have  a  preconceived  notion  of  cost,  or  how  much  effort  he  is 
willing  to  expend  to  arrive  at  a  batter  solution.  Therefore, 
with  soae  knowledge  of  the  grid  optimization  techniques 
available  and  an  understanding  of  the  advantages  and  disad¬ 
vantages  of  each,  the  analyst  can  raalize  the  grid  optimiza¬ 
tion  goal. 

There  are  essentially  two  adaptive  grid  optimization 
strategies: 

(1)  Grid  refineaent,  in  which  the  initial  analysis  is 
performed  on  a  relatively  coarse  grid,  and  new 
degrees-of -freedom  are  added  to  the  same  grid  by  the 
iterative  application  of  the  solution-based 
criterion. 

(2)  Grid  aodif ica tion,  in  which  the  initial  analysis  is 

performed  using  a  prespecified  number  of 

degrees-of-f reedoa,  and  the  solution-based  criterion 
is  employed  to  shift  degrees-of -freedom  from  certain 
regions  to  others.  This  may  involve  complete  grid 
redefinition  in  an  effort  to  obtain  a  near-optimum 
grid  in  a  single  cycle. 


Such  of  the  interest  lately  has  been  in  the  development 
cf  complicated  self-adaptive  software  packages  which  mini- 
size  the  iapact  of  the  user's  skill  on  the  final  solution. 
Ideally,  the  analyst  would  nerely  define  the  problem  and  the 
prograa  would  automatically  generate  and  analyze  the  optimum 
grid  employing  one  or  more  of  these  techniques,  possibly  at 
the  user's  option. 

1.  HATHBHATIC1L  PHOGHAHfllHG 

Ho  discussion  of  grid  optimization  techniques  would  be 
complete  without  a  brief  description  of  mathematical 
programming,  not  only  because  it  is  how  grid  optimization 
was  earliest  attempted,  but  more  importantly,  it  is 
precisely  what  the  engineer  envisions  when  he  hears  the  term 
"optimization''.  It  is  not  a  grid  optimization  technique,  per 
se,  but  rather  a  numerical  process  of  achieving  any  optimi¬ 
zation  objective,  once  it  is  explicitly  defined  in  mathemat¬ 
ical  terms.  In  solid  mechanics  the  finite  element  method  is 
a  numerical  method  for  minimizing  the  potential  energy  func¬ 
tional,  which  in  discretized  fora  may  be  written: 

7T*%uTKu-u^f  (Egn.  4.1) 

where:  £  is  the  global  displacement  vector 

K  is  the  global  stiffness  matrix,  and 
% 

f  is  the  global  load  vector 

In  the  classical  finite  element  formulation,  the  potential 
energy  is  minimized  with  respect  to  the  nodal  displacements, 
which  implies  satisfaction  of  the  following  stationary 
conditions: 


(i  *  1  #  2 ,  •  •  • ,  n) 


(Eqn.  4.2) 


where  n  is  the  number  of  degrees-of- freedom, 
the  very  familiar  systea  of  linear  aquations: 


This  leads  to 


u  -  f  3  0 


(Eqn.  4.3) 


However,  since 


and  f  ar e  functions  of  the  nodal  coordi¬ 


nates,  then  the  potential  energy  could  be  minimized  with 
respect  to  tha  nodal  coordinates  as  well.  This  would  require 
satisfaction  of  the  following  additional  stationary 
conditions: 


(J  8  1,2,..., a) 


(Eqn.  4.4) 


where  a  is  the  number  of  nodal  coordinates,  x^.  This  differ¬ 
entiation  leads  to  the  less  familiar  system  of  non-linear 
equations: 


,  j  3^' 

\  U  ^  n  -  t"*4”  U  =  0 
2  %  3xj  3xj  ^ 


(j  s  1,2,... ,m) 


(Eqn.  4.5) 


This,  then,  is  the  mathematical  statement  of  the  grid  opti¬ 
mization  problem  for  the  elastostatic  case.  The  nodal 
displacement  variables  may  be  eliminated  by  minimizing  the 
potential  energy  with  respect  to  the  nodal  coordinates  only, 
subject  to  the  implicit  constraint  that  Equation  4.3  is 
always  satisfied  [Bef.  4].  Unfortunately  this  does  not  help 
much  because  the  objective  function  is  still  nonlinear, 
rendering  most  numerical  optimization  algorithms  inefficient 
and  unreliable.  The  difficulty  is  even  further  compounded  by 
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the  requirement  that  the  nodal  variables  be  subject  to  side 
constraints  in  order  to  aaintain  the  defined  boundary  of  the 
domain  and  to  ensure  that  the  elenents  neither  distort 
excessively  nor  overlap  one  another.  For  all  except  the 
simplest  of  problems,  these  constraints  may  be  even  more 
severely  nonlinear  than  the  objective  function,  resulting  in 
the  analysis  becoming  prohibitively  expensive  [Ref.  2].  For 
this  reason,  mathematical  programming  in  finite  element  grid 
optimization  has  been  abandoned  in  favor  of  some  equally 
reliable,  yet  far  mere  computationally  efficient  grid  opti¬ 
mization  techniques.  However,  these  early  efforts  with 
mathematical  programming  were  not  totally  in  vain  because 
they  gave  rise  to  the  contouring  techniques. 

B.  COITOORIVS 

Since  mathematical  programming  is  infeasible  for  grid 
optimization,  further  investigations  were  conducted  to 
suggest  some  guidelines  to  enable  the  analyst  to  construct  a 
grid  with  similar  topological  features  of  the  numerically 
optimized  grid  without  the  computational  effort  involved. 
Turcke  [Ref.  4  ],  in  employing  mathematical  programming  in 
the  solution  of  some  simple  two-dimensional  elastcstatic 
problems,  observed  that  there  was  a  very  definite  element 
pattern  commoa  among  problems  involving  high  strain  gradi¬ 
ents  and  that  the  nodes  of  the  numerically  optimized  grid 
generally  tended  to  be  aligned  along  contours  of  some 
response  function  being  modeled.  Consequently,  in  performing 
analyses  on  grids  whose  construction  was  based  on  contours 
derived  from  an  initial  analysis,  it  was  determined  that  the 
following  provided  grid  characteristics  in  regions  of  high 
strain  gradients  similar  to  the  numerically  optimized  grid, 
but  at  a  fraction  of  the  computational  expense: 


•  contours  of  displacement 

•  contours  of  aazieue  principal  stress 

•  contours  of  maximum  shear  stress 

•  contours  of  strain  energy  density  (isoenergetics) 

•  principal  stress  trajectories  (isostatics) 

Since  the  strain  energy  density  is  the  response  which  is 
consistent  with  the  principle  of  ainieua  potential  energy, 
isoenergetics  are  the  most  commonly  used  contours  along 
which  element  edges  are  aligned  [Ref.  4].  However,  there 
still  remains  the  question  of  how  to  position  the  nodes 
along  the  contours.  For  this  reason,  isostatics  have  become 
increasingly  popular  because  the  principal  stress  trajecto¬ 
ries  fora  a  "flow  net"  of  orthogonal  curves  which  can  guide 
the  analyst  in  the  layout  of  the  elements  [Ref.  16]. 

Since  contouring  involves  the  redefinition  of  the  grid, 
as  opposed  to  a  grid  refinement,  its  primary  advantage  is 
that  the  enriched  mesh  is  not  constrained  to  the  element 
configuration  of  the  previous  mesh.  Therefore,  there  is  no 
limit  to  the  amount  of  enrichment  per  cycle  which  can  be 
performed  and  it  is  conceivable  that  an  optimum  mesh  could 
be  generated  in  a  single  cycle  [Bef.  2],  However,  while  the 
computational  cost  of  repeated  analyses  is  reduced,  the 
preprocessing  costs  involved  in  constructing  the  contours 
and  redefining  the  mesh  can  become  quite  high,  especially  if 
the  contours  are  complex.  algorithms  for  performing  these 
tasks  in  two-dimensional  domains  have  been  proposed 
[Bef.  4,  11 ],  but  they  are  not  extendable  to  three- 
dimensional  problems.  The  major  obstacle  for  two-  and 
three-dimensional  domains  is  that  it  is  often  difficult  to 
constrain  the  element  edges  to  th9  contours  without  the 
elements  becoming  elongated  or  distortad  to  the  degree  that 
numerical  inaccuracies  result,  another  difficulty,  not 
addressed  in  the  literature,  is  how  the  contour  increments 
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should  be  selected  when  the  response  function  is  non¬ 
monotonic  over  the  domain. 


C.  SELECTIVE  HEFXIEHEOT 

The  most  commonly  employed  grid  optimization  technique 
is  that  of  selective  refinement.  As  its  name  implies, 
selected  elements  from  a  given  mesh  are  enriched  while  the 
original  element  grid  remains  essentially  intact.  The 
elements  selected  for  refinement  ars  determined  by  the  iter¬ 
ative  application  of  the  solution-based  criterion  to  indi¬ 
cate  which  elements  contribute  most  to  the  solution  error. 
The  refinement  can  be  performed  by  either  the  h-version  or 
the  p-version,  or  even  the  combined  version  if  so  desired, 
but  the  choice  is  most  often  predetermined  by  the  capabili¬ 
ties  of  the  available  preprocessor.  Since  the  addition  of 
new  degrees-of -freedom  over  several  iterations  can  quickly 
enlarge  the  problem,  it  is  advisable  to  perform  the  initial 
analysis  with  a  reasonably  coarse  grid  of  optimally  shaped 
elements,  that  is  nearly  square  quadrilaterals  or  nearly 
equilateral  triangles.  This  is  especially  important  in  the 
h-version  where  it  is  desirable  to  prevent  the  successive 
subdivision  of  elements  from  producing  elongated  new 
elements.  One  refinement  technique  which  will  ensure  this  is 
the  so  called  Mfather-to-f  our  sons"  subdivision  scheme  in 
which  a  single  quadrilateral  or  triangular  element  is 
replaced  by  four  new  ones  by  adding  and  connecting  midside 
nodes  on  the  edges  of  the  original  element  as  shewn  in 
Figure  4.1.  The  major  difficulty  in  selective  refinement 
arises  when  the  addition  of  a  node  along  an  edge  of  the 
element  to  be  subdivided  creates  a  higher  polynomial  ordered 
edge  for  an  adjacent  element  which  is  not  to  be  subdivided. 
There  results  an  incompatibility  in  the  interpolation  of  the 
dependent  variable  along  this  interel9ient  boundary.  Such  is 
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Figure  4.1  Scae  h-Version  Subdivision  Schemes. 

the  case  in  the  h-version  scheme  of  Figure  4.1  and  it  also 
arises  in  the  p- version  when  two  elements  of  different  poly¬ 
nomial  order  share  a  common  edge,  when  this  situation 
occurs,  the  additional  degrees-of- freedom  do  not  actually 
represent  degrees-of-freedom  at  all  because  they  must  be 
numerically  constrained  to  the  polynomial  interpolant  of  the 
lower  order.  Such  constraints  are  usually  imposed  in  one  of 
three  ways:  global  constraint  equations  may  be  written;  the 
constraints  may  be  incorporated  in  the  elemental  basis;  or 
hierarchical  forms  may  be  used  with  the  excess  degrees-of- 
freedom  simply  set  to  zero  and  condensed  out  (Ref.  2]. 

There  are  some  other  selective  refinement  techniques 
which  do  not  require  any  major  software  modifications.  In 
the  h-version,  the  continuity  problem  may  be  circumvented  by 
employing  any  of  the  coarse-to-fina  mesh  transition  schemes 
for  which  all  of  the  element  edges  remain  of  the  same  poly¬ 
nomial  order  [Bef.  17:  p.  210].  However,  it  is  impossible 
to  employ  these  schemes  without  permitting  some  element 
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distortion,  and  the  refinement  must  nearly  always  be 
performed  interactively  rather  than  automatically.  For  the 
inversion,  interelement  continuity  can  be  easily  ensured  by 
employing  var  i  able- node  d  isoparametric  elements,  which 
permit  a  single  element  to  possess  edges  of  differing  poly¬ 
nomial  orders  [Bef.  17:  p.  125}. 

The  analyst  must  also  exercise  care  when  adding  new 
nodes  to  the  boundary  of  the  domain  to  ensure  that  the 
appropriate  boundary  conditions  are  determined  and  applied. 
Furthermore,  if  the  boundary  is  curved,  the  coordinates  of 
the  new  node  should  be  computed  such  that  it  is  placed  on 
the  actual  boundary  and  not  necessarily  on  the  edge  of  the 
element  being  refined  [Bef.  2]. 

The  important  advantage  of  the  selective  refinement 
technigue  is  that  once  an  appropriate  refinement  criterion 
has  been  determined,  selecting  candidate  elements  for 
refinement  in  each  cycle  becomes  straightforward.  The 
refinement  can  then  be  continued  indefinitely  to  achieve 
very  high  accuracy,  but  because  the  solution  phase  is 
repeated  for  each  cycle,  it  is  desirable  to  hold  the  number 
of  cycles  tc  a  minimum.  Because  the  nodes  from  the  previous 
mesh  remain  fixed  for  each  cycle,  selective  refinement  is 
ideally  suited  for  iterative  solution  methods.  The  solution 
values  obtained  from  the  previous  cycle,  combined  with 
interpolated  values  for  the  new  degrees-of-freedom,  provide 
an  excellent  initial  guess  for  the  next  cycle  [Bef.  2]. 

The  major  disadvantage  is  that  the  limited  amount  of 
refinement  which  can  be  performed  in  each  cycle  may  necessi¬ 
tate  several  cycles  tc  obtain  an  optimum  grid.  In  addition, 
if  new  degrees-of-freedom  require  interelement  continuity 
constraints,  data  management  can  become  cumbersome  unless 
the  constraint  is  performed  hierarchically. 


0.  SOBOOBAIH  ISOLATION 


One  of  the  obvious  disadvantages  of  the  selective 
refinement  technique  is  that  the  solution  oust  be  completely 
repeated  for  each  cycle  when,  in  fact,  -he  number  of 
degrees-cf- freedom  added  in  each  cycle  may  be  few  in  compar¬ 
ison  to  the  total  for  the  problem.  In  addition,  the  number 
of  elements  requiring  refinement  in  each  cycle  may  only 
account  for  a  small  portion  of  the  domain.  Although  the 
refinement  criterion  has  indicated  where  the  grid  is  inade¬ 
quate  and  the  approximation  is  likely  to  be  poor,  the  solu¬ 
tion  is  repeated  in  each  cycle  for  those  nodes  where  the 
error  is  presumably  small.  Besides  the  apparent  computa¬ 
tional  inefficiencies,  this  shortcoming  severely  restricts 
the  amount  of  refinement  which  can  be  performed  in  the 
subregions  of  interest  since  it  is  desirable  to  confine  the 
size  of  the  problem  within  reasonable  limits.  An  alternative 
approach  is  to  reformulate  the  problem  for  those  subregions 
where  refinement  is  necessary  and  to  accept  the  results  of 
the  initial  analysis  as  an  adequate  solution  for  the 
remainder  of  the  domain.  The  elements  requiring  further 
refinement,  which  constitute  isolated  subdomains  of  the 
original  problem,  can  generally  be  subjected  to  signifi¬ 
cantly  greater  refinement  than  would  otherwise  be  practical. 
The  solution  obtained  from  the  initial  analysis  is  then  used 
in  imposing  boundary  values  on  those  degrees -of-freedom 
located  on  the  boundaries  of  an  individual  subdomain.  These 
can,  in  turn,  be  used  to  generate  the  boundary  conditions 
for  any  additional  boundary  degrees- of- freedom  introduced  by 
the  refinement  using  an  appropriate  interpolation  scheme. 

This  grid  optimization  technique,  which  the  author  terms 
"subdomain  isolation",  has  some  further  advantages  over 
selective  refinement.  The  subdomain  may  be  selected  arbi¬ 
trarily  small  such  that  excellent  results  may  be  obtained 
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with  a  single  cycle  using  unifora  refinement.  Therefore,  the 
difficulties  involving  coa rse-to-f ine  transition  schemes, 
aleaent  elongation  and  intereleaent  continuity  can  be 
avoided.  Furthermore,  one  can  choose  as  aany  subregions  for 
refineaent  as  desired  without  creating  an  excessively  large 
problem. 

The  obstacle  which  may  prevent  this  technique  from  being 
readily  accepted  is  the  notion  that,  by  imposing  erroneous 
boundary  conditions  on  the  subdomains,  the  convergence  of 
the  finite  element  aethod  to  the  exact  solution  in  these 
regions  has  somehow  been  taapered  with.  This  aversion  may  be 
somewhat  abated  by  considering  a  simple  extension  of 
Saint-Venant *s  Principle  [Ref.  18:  p.  33].  Although  the 

conditions  are  not  rigorously  satisfied  an  the  boundary, 
which  may  result  in  significant  changes  in  the  response 
locally,  the  affect  at  some  sufficient  distance  away  will  be 
negligible.  The  nuaerical  evidence  supports  this  premise. 
While  errors  in  the  boundary  values  may  somewhat  restrict 
the  accuracy  of  the  dependent  variable,  great  improvements 
can  be  realized  in  the  accuracy  of  its  gradients,  which  is 
more  often  the  goal  of  the  optimization.  Since  the  initial 
analysis  provides  the  boundary  values  for  the  subdomains,  it 
is  desirable  that  its  solution  be  as  accurate  as  reasonably 
possible.  Fortunately,  since  subsequent  refinements  are  not 
performed  on  the  original  grid,  the  initial  analysis  may 
involve  significantly  more  degrees-of-freedom  than  in  the 
case  of  selective  refineaent. 

E .  BESH  GRADING 

The  final  grid  optimization  technique  to  be  discussed 
employs  a  mesh  for  which  the  element  sizes  are  successively 
reduced,  according  to  some  geometric  sequence,  towards  a 
selected  region  of  the  domain.  Oae  might  argue  that  mesh 
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grading  is  not  really  an  optimization  technique  since  it  is 
■ost  often  applied  on  an  "a-priori"  basis  rather  than  adap¬ 
tively,  and  that  it  does  not  lend  itself  well  to  the  itera¬ 
tive  application  of  a  solution-based  criterion.  However,  the 
technique  is  siaple  to  use  and  its  ioplenentation  requires 
few  software  aodif ications.  Furthermore,  a  solution-based 
refineaent  criterion  can  be  used  to  give  a  measure  of  the 
quality  of  the  mesh  to  indicate  whether  a  more  pronounced 
grading  aay  prove  beneficial.  Depending  on  the  solution 
parameter  of  interest,  mesh  grading  can  provide  excellent 
accuracy  at  a  low  analysis  cost.  This  refinement  method  must 
therefore  be  considered  among  the  grid  optimization 
techniques. 

For  the  less  elaborate  finite  element  preprocessors, 
aesh  grading  is  often  the  only  refineaent  means  available 
without  resorting  to  a  uniformly  finer  mesh  involving  many 
■ore  degrees-of-freedca.  The  aost  common  method  of  implemen¬ 
tation  in  two-diaensions  is  to  first  define  the  problem 
domain  in  tens  of  a  curvilinear  quadrilateral  by  selecting 
four  keynodes  along  the  problea  boundary.  Then  the  boundary 
nodes  are  spaced  according  to  soae  geometric  sequence  based 
on  the  user-provided  bias  parameters  which  determine  the 
density  of  the  nodes  towards  selected  points  on  the  four 
quadrilateral  edges.  Finally,  curves  are  generated  to 
connect  the  boundary  nodes  on  opposite  edges  of  the  quadri¬ 
lateral,  thus  producing  a  graded  aesh.  This  process,  which 
can  also  be  extended  to  three-dimensions,  is  the  mesh  gener¬ 
ation  scheae  employed  in  the  finite  element  code  GIFTS 
[fief.  19]. 

The  major  disadvantage  of  aesh  grading  is  that  in  order 
to  achieve  sufficiently  small  elaments  in  the  region  of 
interest,  the  elements  must  grow  successively  larger  away 
froa  that  region.  This  aay  be  very  undesirable,  especially 
if  refineaent  is  called  for  in  aora  than  one  region  of  the 
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doaain,  in  which  case  the  aesh  aust  be  generated  and  graded 
by  subdoaains,  thereby  coaplicating  the  data  management 
involved. 

Another  disadvantage  is  that  unless  the  doaain  possesses 
soae  special  geometric  symmetry,  excessive  element 


Figure  *.2  Graded  Hesh  for  a  Perforated  Square  Plate. 

elongation  will  usually  result  if  a  highly  pronounced 
grading  is  required.  Soae  configurations  are  particularly 
well  suited  for  refinement  for  aesh  grading  such  as  the 
classical  perforated  square  plate  problem  in  solid  mechanics 
shown  in  Figure  4.2. 
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Now  that  the  necessary  tools  for  performing  grid  optimi¬ 
zation  have  been  introduced ,  it  is  time  to  employ  them  in  an 
attempt  to  obtain  optimal  solutions  to  some  practical  prob¬ 
lems  in  engineering.  An  obvious  starting  point  for  such  an 
investigation  is  the  one-dimensional  boundary  value  problem. 
Bhile  most  of  the  fruitful  research  in  grid  optimization  has 
concentrated  on  problems  of  higher  dimensions,  the  one- 
dimensional  problem  is  a  very  convenient  device  for  studying 
finite  element  grid  optimization.  Foremost,  one-dimensional 
finite  element  models  possess  a  unique  connectivity  in  that 
adjacent  elements  meet  at  their  end  nodal  points.  Therefore, 
refinement  by  the  h-  or  p-versions,  or  by  relocating  nodal 
points  becomes  a  trivial  task,  which  does  not  involve  any  of 
the  difficulties  so  frequently  encountered  with  higher 
dimensional  problems,  such  as  preserving  interelement  conti¬ 
nuity  and  maintaining  optimal  element  shapes.  Furthermore, 
one- dimensional  studies  can  often  provide  valuable  insight 
to  the  solution  of  more  difficult  higher  dimensional 
problems. 

The  primary  concerns  in  the  selection  of  the  problems  to 
be  studied  were  as  fellows: 

(1)  there  should  exist  an  analytical  solution  to  provide 
a  means  of  reliable  error  analysis; 

(2)  the  solution  and  its  resultants  should  exhibit 
sufficiently  high  gradients  so  that  the  effective¬ 
ness  of  the  grid  optimization  is  readily  observed. 
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Because  of  the  complexity  and  a  certain  degree  cf  arbi¬ 
trariness  involved  in  the  conputation  of  element  residuals 
and  a-posteriori  error  estimates,  the  solution  parameter 
variation  is  the  refinement  criterion  of  choice.  There  are 
several  solution  parameters  which  are  easily  computed, 
requiring  minimal  software  changes  to  an  existing  finite 
element  code. 

Furthermore,  for  the  one-dimensional  investigation,  it 
was  decided  to  simplify  the  analysis  by  exploiting  the 
linear  elements,  ihile  it  is  granted  that  improved  solution 
accuracy  may  generally  be  obtained  by  employing  higher  order 
elements,  it  will  be  assumed  that  conclusions  based  on  the 
use  of  linear  elements  can  be  applied  as  well  to  elements  of 
higher  polynomial  order. 


A.  ELASTIC  CABLE  PROBLEB 

Consider  an  elastic  cable  under  tension  r,  stretched 
between  two  points  a  distance  2 L  apart.  If  the  cable  is 
supported  by  a  Rinkler,  or  elastic,  foundation  of  modulus  k, 
and  a  concentrated  load  P  is  applied  at  the  midspan,  the 
resulting  deflection  v(x)  ,  (0  <  x  <  L)  ,  is  as  shown  in 

Figure  5.1.  The  analytical  solution  and  finite  element 
approach  for  this  problem  are  presented  in  Appendix  A. 

The  initial  finite  element  analysis  was  performed  using 
ten  linear  elements  of  uniform  length.  From  this  initial 
analysis  the  approximate  distribution  of  the  following  solu¬ 
tion  resultants  was  obtained  over  the  domain: 

•  the  displacement,  v  (the  solution) 

•  the  slope,  v* 

•  the  strain  energy,  0 

•  the  strain  energy  density,  SED  (du/dx) 


Figure  5.1  Tension  Cable  Deflection  on  Elastic  Foundation 


Subsequent  analyses  were  performed  for  finite  element  models 
using  the  same  number  of  elements,  but  with  the  nodes  redis¬ 
tributed  to  achieve  approximately  uniform  variation  of  the 
above  parameters  over  each  element.  Note  that  the  strain 
energy  refinement  criterion  producas  elements  of  identical 
strain  energy  content.  In  addition,  the  problem  was  solved 
employing  graded  element  models  of  various  adjacent  element 
length  ratios.  The  resulting  element  models  based  on  these 
refinement  criteria  are  shown  in  Figure  5.2  (a-f).  The 
graded  model  (b)  for  an  element  length  ratio  of  1.2  is 
presented  for  comparison  because  it  produced  good  overall 
solution  results. 

is  previously  mentioned,  the  solution  resultants  of 
primary  interest  are  the  maximum  displacement,  the  maximum 
slope,  and  the  integral  of  the  displacement  over  the  domain, 
because  they  represent  important  analogous  solution  results 
in  nearly  all  fields  in  which  finite  element  analysis  is 
often  performed.  The  accuracy  obtained  in  these  values  for 
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Pigure  5.2  Tension  Cable  Refinements. 


4  1 


each  of  the  refined  models  is  presented  in  Table  I.  As  car. 
be  seen  in  Figure  5.2,  the  strain  energy  and  strain  energy 
density  criteria  produced  extreme  variations  of  element 
length  while  the  criteria  of  displacement  and  slope  result 


TABLE  I 

Tension  Cable  Problen  solution  Results 


Problem  Parameters:  L  = 

k  = 
T  = 
P  = 

1 00  in 

1  psi 
1000  lb 
1000  lb 

Variation 

Befinement 

Criterion 

Percentage  Relative 

v  (max)  v'  (max) 

Error 

Uniform 

-0.40 

0.36 

0.12 

Graded  (1.  2) 

-0. 19 

0.  07 

0.17 

V 

-0.  18 

0.  06 

0.39 

V' 

-0.23 

0.05 

0.87 

0 

-1.03 

0.05 

3.58 

SED 

-1.  29 

0.05 

4.  17 

0  (mod) 

-0.  53 

0.04 

1.63 

SED  (mod) 

-0.51 

0.03 

1.48 

in  more  moderate  variations.  It  can  be  observed  in  Table  I 
that  the  more  pronounced  refinements  based  on  energy  distri¬ 
bution  result  in  greater  accuracy  for  the  maximum  slope  but 
with  the  accompanying  severe  penalty  of  significantly  poorer 
estimations  nf  the  maximum  displacement  and  the  integral 
quantity.  For  this  particular  problem  the  uniform  grid 
provides  optimal  accuracy  of  the  integral  quantity. 
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therefore  refinement  cannot  reduce  its  error.  Yet  great 
improvement  in  the  accuracy  of  the  maximum  slope  and  modest 
improvement  in  the  accuracy  of  the  maximum  displacement  can 
be  achieved  with  moderate  refinements  based  on  the  displace¬ 
ment  and  slope  distributions. 

One  might  assume,  and  correctly  so,  that  the  ability  of 
the  energy  refinements  to  produce  the  best  accuracy  for  the 
maximum  slope  is  due  to  the  extremal/  small  elements  which 
result  in  the  area  where  that  quantity  occurs.  Furthermore, 
it  would  be  correct  to  propose  that  the  reason  for  these 
refinements  producing  poorer  estimates  than  the  uniform 
model  for  the  other  two  quantities  of  interest  is  that  the 
excessively  large  elements  in  the  ragions  of  low  gradients 
severely  overstiffen  the  model  there.  It  would  then  seem 
plausible  tc  improve  the  accuracy  for  the  maximum  displace¬ 
ment  and  the  integral  quantity  by  redistributing  the  nodes 
in  these  regions  to  prevent  such  excessively  large  elements. 
This  was  done  by  arbitrarily  employing  a  grading  scheme  to 
the  three  largest  elements  to  produce  the  modified  refine¬ 
ments  based  on  strain  energy  and  strain  energy  density  shown 
in  Figure  5.2  (g)  and  (h) .  As  can  be  seen  in  Table  I,  such  a 
modification  did  indeed  significantly  reduce  the  errors  in 
the  maximum  displacement  and  the  integral,  but  it  even 
further  improved  the  accuracy  for  the  maximum  slope. 

One  might  conclude  from  Table  I  that  the  best  overall 
model  was  obtained  using  the  graded  mesh,  and  that  since  it 
is  easier  to  obtain,  it  should  be  deemed  the  optimal  grid. 
But  this  particular  grading  was  caosen  for  precisely  that 
reason  and  was  presented  only  as  a  means  of  comparison.  In 
practice,  the  selection  of  a  grading  ratio  is  somewhat  arbi¬ 
trary  and  making  an  adequate  choice  may  be  difficult. 

There  is  justifiable  confusion  as  to  which  refinement 
produced  the  '•best'*  solution  accuracy  for  this  problem  and 
it  raises  perhaps  the  most  important  issue  in  the  subject  of 
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grid  refinement.  Before  any  optimization  process  can  be 
pursued,  the  optiaization  goals  oust  be  explicitly  defined. 
Clearly,  as  is  the  case  in  this  problem,  the  designation  of 
the  optimum  grid  would  depend  heavily  upon  which  of  the 
three  solution  resultants  is  most  critical  to  the  analysis. 

B.  TAPEBBD  BAB  PBOBLB9 

The  linearly  tapered  bar  under  axial  loading  has 
received  considerable  attention  and  was  one  of  the  early 
problems  for  which  analytical  grid  optimization  was 
employed.  Consider  a  tapered  elastic  bar  of  length  L  and 
modulus  E,  fixed  at  one  end,  with  an  axial  load  P  applied  at 
the  other,  .  for  which  the  axial  displacement  u(x), 
(0  <  x  <  L)  ,  is  desired.  The  cross-sectional  area  varies 
linearly  from  A  at  the  fixed  end  to  At  at  the  tip,  as  shown 


Figure  5.3  Linearly  Tapered  Bar  Under  Axial  Loading. 


in  Figure  5.3.  The  analytical  solution  and  finite  elemen 
approach  are  presented  in  Appendix  B. 


One  of  the  significant  features  of  the  tapered  bar 
problea  is  that  the  aaziaua  stress  can  be  very  difficult  to 
aodei  accurately,  and  it  is  for  precisely  these  problems 
exhibiting  large  strain  gradients  that  grid  optimization 
becoaes  aost  beneficial.  Interestingly,  the  stresses 
obtained  at  the  eleaent  midpoints  are  exact  f  '■  *  this 
problea,  and  the  difficulty  arises  froa  the  inability  of  the 
constant  slope  shape  functions  to  model  the  maxiaum  stress 
occurring  at  the  boundary.  In  exaaining  this  problem, 
Prager  [Hef.  20]  deaonstrated  analytically  that  when  each 
eleaent  has  the  sane  strain  energy  content,  the  relative 
error  in  displaceaent  is  identical  for  all  the  nodes. 
However,  this  phenoaenon  appears  peculiar  to  this  problem 
and  the  author  does  not  subscribe  to  such  a  measure  of  an 
optimum  grid.  Judging  the  effectiveness  of  a  particular 
refineaent  based  upon  the  deviation  or  the  Bean  value  of  the 
pointvise  errors  generally  tends  to  be  unfavorable  to  opti¬ 
mization  procedures  since  they  almost  always  introduce  many 
aora  nodes  in  those  regions  where  the  response  is  aost 
difficult  tc  aodei.  Hence,  an  improved  solution  may  have  a 
larger  aean  value  of  the  pointwise  errors  [Kef.  3]. 

The  criteria  eaployed  in  the  refineaent  of  the  tapered 
bar  aodei  are  identical  to  those  used  in  the  cable  problem 
and  their  effects  are  shown  in  Figure  5.4  (a-e) .  Two  excep¬ 
tions  are  that  now  the  displaceaent  and  strain  energy 
criteria  produce  identical  refineaents,  and  the  graded  model 
chosen  as  the  best  overall  is  now  based  on  a  grading  ratio 
of  1.4,  producing  a  more  drastic  refinement  than  that  of  1.2 
for  the  cable.  This  further  demonstrates  the  difficulty 
involved  in  obtaining  adequate  eleaent  grids  on  an 
"a-priori"  basis. 

The  solution  results  are  presented  in  Table  II  and  the 
aost  readily  apparent  observation  for  this  problem  is  the 
large  errors  in  the  maxiaum  slope,  which  would  severely 
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(b)  Graded  1.4 


TABLE  II 

Tapered  Bar  Problee  Solution  Basalts 


Problem  Pazaaeters: 


L 

A 

A 

E 

P 


3?? 


13 
0.5 
10x10* 
10x103 


£n 

in2 

in* 

IS1 


Variation 

Refineaent 

Criterion 

Percentage  Relative 

u(max)  u*  (aax) 

Error 

^udx 

Onifora 

-3.80  -37.5 

0.68 

Graded  (1.  4) 

-0.78  -4.1 

0.14 

u  ;  U 

-0.85  -10.6 

0.15 

u* 

-1.81  -7.7 

0.33 

SED 

-6.54  -3.6 

1.18 

SED  (aod) 

-1.99  -3.6 

0.36 

onderestieate  the  naxlaaa  stress.  These  results  are  based  on 
quadratic  extrapolation  of  the  exact  slopes  at  the  element 
midpoints,  since  the  linear  shape  functions  would  produce 
even  poorer  estiaatiens  of  the  aaxiauo  slope.  As  before,  the 
aore  extreae  refineaent  based  on  the  strain  energy  dansity 
variation  provides  the  aost  accurate  estiaation  of  the 
aaxiaua  slope,  but  with  the  accompanying  degradation  in 
estiaates  for  aaxiaua  displacement  and  the  integral  of 
displaceaent.  Again,  the  large  errors  in  these  values  may 
be  significantly  reduced  by  employing  a  grading  scheme  to 
restrict  the  size  of  the  larger  elements  as  shown  in  Pigure 
5.4  (f)  .  Unlike  the  previous  problem,  such  a  modification 
has  no  effect  on  the  estimate  of  maximum  slope  because  of 
the  extrapolation  of  the  element  midpoint  slopes,  which  are 
exact  regardless  of  the  element  model. 


&  different  version  of  the  tapered  bar  problem,  tor 
which  the  displacement  and  strain  energy  criteria  will  not 
produce  identical  refinements,  involves  replacing  the 
concentrated  tip  load  P  with  a  linearly  varying  axially 
distributed  load  g(x),  specified  by  the  values  at  the  fixed 
end  ^  and  the  tip  gt .  The  problem  may  be  further  modified 
by  reversing  the  bar  such  that  the  maximum  slope  occurs  at 
the  fixed  end,  while  the  maximum  displacement  occurs  at  the 


Figure  5.5  lever sed  Tapered  Bar  with  Distributed  Load. 


free  end  as  shown  in  Figure  5.5.  The  case  of  the  linearly 
varying  distributed  load  is  included  in  the  formulation  in 
Appendix  B. 

This  problem  was  solved  for  a  uniformly  distributed  load 
using  the  same  procedure  as  in  the  previous  two  problems. 
The  refinement  models  are  presented  in  Figure  5.6  and  the 
solution  results  in  Table  III.  Tha  observations  are  consis¬ 


tent  with  those  made  in  the 
would  likely  agree  that  th 


previous  problems,  but  now  one 
e  refinement  based  on  the  strain 


TABLE  III 

Reversed  Tapered  Bar  solution  Results 


Problem  Parameters: 


100 

0.5 

10.5 

10x10* 

100 


psi 

lb/ii 


Variation 

Refinement 

Criterion 

Percentage  Relative 

u  (max)  u1  (max) 

Error 

Sa; 

Uniform 

-5.5  -39.4 

-7.1 

u 

-2.0  -7.5 

-2.6 

u' 

-2.7  -8.1 

-3.4 

u 

-3.7  -5.9 

-4.7 

SED 

-11.9  -3.4 

-15.3 

SED  (mod) 

-3.1  -3.4 

-4.0 

energy  density  would  represent  an  optimal  grid,  provided 
that  modifications  are  introduced  to  prevent  any  elements 
from  growing  excessively  large. 


C.  GUIDELINES  FOB  ONE-DIMENSIONAL  GRID  OPTIMIZATION 

The  most  important  lesson  to  be  learned  from  this  one- 
dimensional  study  is  that  the  grid  optimization  procedure  is 
necessarily  dictated  by  the  optimization  goal,  or  the  under- 
lying  purpose  for  performing  the  finite  element  analysis.  No 
element  grid  can  possibly  provide  optimum  accuracy  for  every 
solution  resultant  cf  interest.  In  solving  these  simple 
problems,  a  balance  has  been  sought  for  achieving  adequate 
accuracy  for  three  of  the  more  important  solution 


resultants,  with  emphasis  on  the  maximum  value  of  the  deriv¬ 
ative  of  the  dependent  variable,  which  acre  often  is  not 
only  the  acst  important  part  of  the  solution  but  also  the 
most  difficult  to  obtain  accurately  in  finite  element 

analysis. 

The  important  grid  optimization  techniques  of  intro¬ 
ducing  mere  degrees-of-freedom  by  subdividing  the  elements 
or  increasing  their  polynomial  order  have  been  intentionally 
omitted  in  favor  of  the  optimization  strategy  of  seeking 
maximum  solution  accuracy  for  a  specified  number  of 
degrees-of-£reedom  using  linear  elements.  This  is  because 
such  a  procedure  is  not  so  straightforward  in  two- 
dimensional  problems  where  the  number  of  degrees-of-freedom 
are  dependent  on  some  geometric  considerations,  which  do  not 
appear  in  problems  of  one-dimension.  Based  on  this  choice  of 
optimization  strategy,  it  appears  the  strain  energy  density 
variation  provides  the  most  useful  criterion  for  the  adap¬ 
tive  refinement  of  the  initial  grid,  ret  all  three  problems 
demonstrated  some  pathological  results  that  can  arise  when 
the  elements  are  permitted  to  grow  excessively  large  in  the 
regions  where  the  strain  energy  density  varies  the  least,  in 
applying  a  scheme  to  restrict  the  size  of  the  largest 
elements,  no  mention  has  been  made  of  how  tc  determine  when 
an  element  is  excessively  large.  It  has  become  the  experi¬ 
ence  of  the  author  that  any  element  representing  over  half 
of  the  domain  should  probably  be  considered  too  large,  and 
measures  should  be  employed  to  restrict  its  size. 

It  would  appear,  at  least  for  these  classes  of  problems, 
that  this  difficulty  of  decreasing  accuracy  of  a  particular 
solution  parameter  for  successive  refinements  can  be  ignored 
by  merely  accepting  the  largest  value  among  the  cycles  as 
the  most  accurate  solution  result.  For  example,  it  was 
demonstrated  that  the  refinement  based  on  strain  energy 
density  provided  significant  improvement  in  the  accuracy  for 
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the  laziaua  3lope  but  underestimated  the  maximum  displace¬ 
ment  even  more  than  the  initial  uniform  grid.  Assuming  that 
the  linear  element  model  always  underestimates  such  maxima, 
the  maximum  slope  for  refined  grid  and  the  maximum  displace¬ 
ment  for  the  unrefined  grid  could  be  accepted  as  the  optimal 
results  of  the  analysis.  The  fallacies  of  such  a  procedure 
are  that,  first,  the  refinement  may  not  represent  the 
optimal  grid  as  it  has  been  defined  and,  second,  for  self- 
adaptive  finite  element  codes  the  user  is  provided  with  the 
"optimum  grid"  of  the  final  cycle  and  the  solution  results 
thereof. 

Helosh  and  Harcal  [Bef.  21]  have  proposed  an  alternative 
use  of  the  refinement  criterion  based  on  strain  energy 
density  variation  which  avoids  the  problem  of  excessive 
element  growth  altogether.  Beginning  with  a  reasonably 
coarse  unifori  grid,  those  elements  with  the  greatest  strain 
energy  density  variation  are  selectively  refined  by  either 
subdividing  them  or  increasing  their  polynomial  order  with 
the  introduction  of  additional  degrees-of-freedom.  while 
such  a  procedure  does  not  equi-distribute  the  element  strain 
energy  variations,  it  can  reduce  them  all  to  some  prespeci¬ 
fied  tolerance,  such  as  a  percentage  of  the  average  element 
variations  for  the  initial  analysis.  Because  this  procedure 
is  particularly  attractive  for  gril  refinement  in  problems 
of  higher  dimensions,  it  will  be  employed  extensively  for 
the  study  of  grid  optimization  for  two-dimensional  problems 
in  the  next  chapter. 
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Since  investigators  began  working  in  the  field  of  finite 
element  grid  optimization  in  the  early  1970's,  nearly  all  of 
the  effort  has  been  devoted  to  the  development  of  a  system¬ 
atic  procedure  for  obtaining  optimal  grids  for  two- 
dimensional  problems  of  elasticity.  Even  today  there  are 
several  competing  approaches  to  this  problem  and  no  partic¬ 
ular  one  has  yet  been  overwhelmingly  accepted  as  the 
preferred  method  of  grid  optimization.  While  it  is  the  two- 
dimensional  problem  for  which  most  of  these  techniques  have 
been  developed,  their  application  to  such  can  be  much  more 
difficult  than  for  the  one-dimensional  case.  Almost,  invari¬ 
ably  when  performing  grid  refinement  on  two-dimensional 
domains,  the  analyst  is  confronted  with  the  problems  of 
maintaining  interalement  compatibility  and  preventing  severe 
element  distortion. 

In  selecting  an  appropriate  two-dimensional  problem  for 
the  application  of  some  grid  optimization  techniques  and  a 
comparison  of  their  effectiveness,  it  is  desirable  that  the 
test  case  possess  the  following  properties: 

•  the  analytical  solution  should  exist  in  order  to 
perform  reliable  error  analysis 

•  the  solution  should  exhibit  sufficiently  large 
gradients  to  provide  a  meaningful  measure  of  the 
refinement  effectiveness 

•  the  idealization  should  have  one  degree-of-freedom  per 
node  and  possess  simple  boundary  conditions  to 
minimize  the  computational  effort  involved  in 
repeated  solutions 
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There  are  few  problems  that  meat  these  criteria,  but 
Saint-?enant  torsion  of  a  non-circular  section  provides  a 
good  test  case. 

A.  PROBLEM  DESCRIPTION 

Consider  a  solid  circular  shaft  of  radius  "a"  made  from 
isotropic  material  of  shear  modulus  s  and  having  a  circular 
groove,  or  keyway,  of  radius  b  along  a  generator  of  the 
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shaft.  The  shaft  cross-section  is  shown  in  Figure  6.1.  The 
shaft  is  subjected  to  an  applied  torque  T  which  produces 
an  angle  of  twist  per  unit  length  9.  The  problem  may  be 

solved  by  finding  the  Prandtl  torsional  stress  function  ^ 
which  satisfies  the  governing  differential  equation: 

.2  +  s,,2  +2  =  0  (Eqn.  6.1) 
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subject  to  the  Dirichlat  condition  that  ♦  *  0  on  the 
section  boundary.  The  torsional  stress  function  is  defined 
such  that  the  shear  stress  t  at  any  point  on  the  domain  may 
be  expressed  as: 


T  -G0  [(  |f  l2  ♦  (  fj  l2]^  (Em-  6.2) 

For  this  formulation,  the  angle  of  twist  9  is  prescribed, 
rather  than  the  applied  torgue  T.  The  torque  is  is  calcu- 
lated  from  the  area  integral: 


(Eqn.  6.3) 


The  analytical  solutions  of  Equations  6.1  and  6.2  are 
derived  by  Sokolnikoff  [Eef.  22:  pp.  1*1- 1*3]  and  are 
presented  in  Appendix  C  along  with  the  evaluation  of 
Equation  6.3  and  a  prescribed  finite  element  formulation. 
For  this  problem,  the  three  solution  resultants  of  interest 
for  the  grid  optimization  study  are: 

(1)  maximum  value  cf  the  dependent  variable,  or  torsion 
function  ^iax; 

(2)  maximum  value  cf  the  gradient  of  the  dependent  vari¬ 
able  (a  quantity  proportional  to  maximum  shearing 
stress  T max)  ; 

(3)  the  area  integral  of  the  dependent  variable  over  the 
domain  (a  quantity  proportional  to  the  applied 
torque  T) . 

These  quantities  -  the  dependent  variable,  its  gradient,  and 
an  integral  thereof  -  are  selected  as  representative  of 
entities  whose  error  one  might  wish  to  minimize  in  a  finite 

element  analysis. 
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is  can  be  seen  in  Figure  6.1 ,  the  donain  of  this  problem 
is  symmetric  about  the  x-axis,  therefore  the  finite  element 
solution  need  only  be  obtained  for  the  upper  half  of  the 
doaain.  For  all  of  the  solutions  presented  herein  the 
probles  geoaetry  is  defined  by  assigning  the  diaensionless 
ratio,  b/a  *3.4,  and  an  acceptable  upper  limit  on  the  anal¬ 
ysis  cost  was  arbitrarily  chosen  to  be  that  corresponding  to 
approxiaately  one  hundred  nodal  points.  The  computation  and 
assembly  of  the  finite  eleaent  matrices  and  solution  of  the 
resulting  systea  of  equations  was  performed  using  the  steady 
state  heat  conduction  operations  of  CAL-NPS  [  Hef .  23].  This 
group  of  subroutines  coaprises  an  efficient  finite  eleaent 
code  for  solving  Poisson's  eguation  in  two  or  three  dimen¬ 
sions  and  has  the  additional  advantage  of  permitting 
variable-noded  isoparametric  eleaents. 

Since  there  was  no  readily  available  interactive  prepro¬ 
cessor  which  lent  itself  well  to  adaptive  mesh  refinement, 
the  author  had  no  choice  but  to  create  his  own.  Since  the 
problea  doaain  is  siaply  connected,  the  automatic  mesh 
generation  was  performed  employing  inverse  mapping  of  a 
single  cubic  isoparametric  eleaent  of  the  serendipity  family 
onto  the  problea  dcaain  [Bef.  24:  pp.  228-229].  Napped 
boundary  nodes  were  repositioned  to  conform  to  the  actual 
doaain  boundary  and  additional  nodes  generated  during  the 
refineaent  process  were  mapped  using  the  same  procedure. 

Since  the  finite  eleaent  code  selected  for  this  investi¬ 
gation  provided  output  only  for  the  nodal  values  of  the 
dependent  variable,  it  was  coupled  to  the  author's  postpro¬ 
cessor.  Such  a  postprocessor  is  necessary  in  the  optimiza¬ 
tion  process  for  computing  nodal  values  of  shear  stresses 
and  strain  energy  density,  eleaent  contributions  to  torque 
and  total  strain  energy,  and  exact  results  from  theory. 


C.  ASYMPTOTIC  EHfiOH  ANALYSIS  FOR  ONIFORH  EEFTHEHENT 


The  concept  of  asymptotic  convergence  rate  for  uniformly 
refined  grids  was  presented  in  Chapter  3.  When  the  number  of 
uniformly  distributed  deg rees-of-f readom  is  sufficiently 
large,  the  log-log  plot  of  th9  relative  energy  error  versus 
the  number  of  degrees-of-f reedom  is  approximately  linear  in 
the  asymptotic  range.  The  slope  of  this  line  represents  the 
asymptotic  rate  of  convergence  in  energy. 

It  so  happens  that  relative  error  in  the  torque  T  of 
this  problem  is  equal  to  the  relative  energy  error  and 
therefore  exhibits  this  linear  asymptotic  behavior  on  the 
log-log  plot  against  the  number  of  uniformly  distributed 
degrees-of- freedom.  Fortunately,  the  other  two  solution 
resultants  of  interest  behave  similarly.  This  will  prove 
very  beneficial  in  performing  the  error  analysis  for  this 
two-dimensional  study  for  two  reasons.  First,  because  it  is 
unnecessarily  difficult  to  construct  an  optimal  grid  with 
the  same  number  of  degrees- of- freedom  as  a  uniform  grid,  the 
linear  behavior  of  the  solution  resultants  in  the  asymptotic 
range  on  the  log-log  scale  permits  interpolation  for  any 
number  of  degrees-of -freedom.  Then  the  solution  results  for 
a  uniform  grid  of  the  identical  number  of  degrees-of-f  reedom 
provides  a  reference  for  comparison  to  determine  the  effec¬ 
tiveness  of  the  optimization  technique.  Secondly,  if  the 
convergence  rate  of  a  particular  solution  resultant  is 
extremely  slow,  as  is  often  the  case  for  maximum  stress,  it 
becomes  difficult  to  gain  an  appreciation  of  the  true  effec¬ 
tiveness  of  the  optimization.  For  example,  an  order  of 
magnitude  reduction  in  the  relative  solution  error  may 
require  an  order  of  magnitude  increase  in  the  number  of 
degrees-of- freedom  using  uniform  refinement,  but  relatively 
few  additional  degrees-of -freedom  using  an  optimization 
technique.  Therefore,  it  will  be  enlightening  to  extrapolate 


the  relative  error  versus  degrees-of- freedom  curve  to  obtain 
a  rough  approximation  of  the  number  of  degrees-of-f reedom 
necessary  to  obtain  solution  accuracies  similar  to  the 
optimal  grid,  but  using  successively  finer  uniform  grids.  Of 
course,  this  is  only  an  estimation  and  ignores  such  reali¬ 
ties  as  numerical  ill-conditioning  and  computer  round-off 
error. 

A  uniform  grid  is  one  for  which  all  of  the  elements  are 
of  the  same  size  h  and  the  same  polynomial  order  p.  Clearly, 
it  is  impossible  to  obtain  such  a  grid  for  this  particular 
domain  using  isoparametric  mapping,  but  a  nearly  uniform 
grid  may  be  constructed  in  which  the  elements  are  of  approx¬ 
imately  the  same  size,  such  uniform  grids  are  shown  for  the 
cases  of  linear  quadrilateral  elements  and  quadratic  seren¬ 
dipity  elements  (Pig.  6.2).  For  this  geometry,  the  uniform 
grid  is  not  uniquely  defined  for  a  specified  number  of 
elements.  This  is  because,  in  performing  isoparametric 
mapping,  there  must  be  specified  four  Iteynodes  on  the  actual 
domain  boundary  to  correspond  to  the  four  corner  nodes  of 
the  parent  square.  Since  this  domain  has  only  three 


Figure  6.3  Keynode  Placement  for  Isoparametric  Sapping. 
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vertices,  the  placement  of  the  fourth  keynode  is  at  the 
discretion  cf  the  analyst  (Fig.  6.3)  ,  and  can  have  a  notice¬ 
able  effect  on  the  solution  results. 

The  asyaptotic  error  analysis  was  performed  for  the 
three  solution  resultants  of  interest  using  uniform  grids  of 
linear  and  quadratic  elements.  The  results  are  presented  in 
Figure  6.4.  ill  of  the  solution  resultants  behaved  as 
predicted  with  the  exception  of  the  maximum  torsion  function 
value  using  linear  elements,  U  (1)  .  It  appears  that  the 
accuracy  of  this  particular  parameter  is  very  strongly 
dependent  upon  the  keynode  placement.  The  curve  constructed 
in  Figure  6.4  represents  an  average  for  several  keynode 
positions. 

Hhile  Figure  6.4  is  intended  primarily  to  serve  as  a 
reference  tool  for  future  analyses,  it  provides  some  inter¬ 
esting  irformation: 

(1)  For  the  cases  of  maximum  torsion  function  value  and 
applied  torque  (and  energy)  ,  the  asymptotic  rate  of 
convergence  using  quadratic  elements  is  more  than 
twice  that  for  linear  elements. 

(2)  Bhile  the  error  in  torque  for  the  quadratic  case  is 
always  smaller  than  that  for  the  linear  case,  the 
linear  grid  may  provide  better  accuracy  for  the 
maximum  torsion  function  value  V^max  in  the  pre- 
asymptotic  range. 

(3)  Both  accuracy  and  convergence  rate  in  the  maximum 
shear  stress  are  only  minutely  greater  for  the  quad¬ 
ratic  element  grid  than  for  the  linear  one. 

However,  for  this  last  observation,  the  point  must  be  made 
that  for  the  linear  element  grid,  th9  maximum  shear  stress 
was  obtained  by  quadratic  interpolation  rather  than  from  the 
linear  shape  functions,  While  this  will  greatly  improve  the 
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accuracy  cf  the  aaxiaua  shear  stress  approximation,  it  will 
have  no  effect  on  its  rate  of  convergence.  Therefore,  if 
obtaining  an  optiaal  estiaate  of  the  aaxiaua  shear  was  the 
purpose  of  the  analysis,  there  is  such  to  be  said  on  behalf 
of  linear  eleaents  besides  their  computational  efficiency. 
Of  course,  this  observation  is  based  on  unifora  grid  refine- 
ment,  which  would  rarely  coapete  favorably  with  the  optiai- 
zation  technigues  to  be  exaained. 

The  reason  that  the  rate  of  convergence  in  maximum  shear 
stress  is  so  poor  using  uniform  refinement  for  this  problem 


TC  A)=0 
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Figure  6.5  Stress  Distribution  on  Shaft  Keyway. 

can  be  seen  in  Figure  6.5.  The  shear  stress  varies  greatly 
over  a  short  distance,  by  increasing  from  zero  at  point  &  to 
its  aaxiaua  value  at  point  B.  Is  a  result,  there  exists  a 
region  of  excessively  large  strain  gradients  along  the 
keyway  which  severely  hinders  the  rate  of  convergence  when  a 
unifora  grid  is  eaployed.  If  the  keyway  radius  were  allowed 
to  approach  zero  producing  a  singularity  in  the  solution,  it 
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would  likely  be  necessary  to  employ  even  higher  order 
elements  via  the  p -vers ion  in  order  to  achieve  convergence 
using  uniform  refinement  [Ref.  1]. 

D.  PROBLEB  SOLOTIOH  1ITH  GRID  OPTIBIZATIOH 

The  finite  element  solution  of  the  torsion  problem  will 
be  obtained  employing  the  following  grid  optimization  tech¬ 
niques  as  presented  in  Chapter  4: 

•  Contouring 

-  contours  of  the  torsion  function;  linear  elements 

-  contours  of  shear  stress;  linear  elements 

-  contours  of  strain  energy  density;  linear  elements 

•  Selective  Refinement 

-  h-version;  linear  elements 

-  h-version;  quadratic  elements 

-  p- version 

•  Subdcmain  Isolation 

-  linear  elements 

-  quadratic  elements 

•  Besh  Grading 

-  linear  elements 

-  quadratic  elements 


* 

The  original  finite  element  analysis  was  performed 
on  a  uniform  grid  of  98  linear  ele«9nts,  78  nodes,  and  72 
degrees-of-fraedom.  The  finite  element  solution  provided  the 
nodal  values  of  the  torsion  function  ,  from  which  the 
conventional  nodal  resultants  of  shear  stress  t  ,  and 
applied  torque  T  were  computed.  Based  upon  the  maximum  and 
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minimum  values  obtained  for  each  parameter,  along  with 
consideration  for  their  values  along  the  boundary,  the 
contours  to  be  used  for  nodal  plasaaent  in  each  case  were 
selected.  The  nuaber  of  contours  for  each  case  was  chosen  to 
aaintain  approximately  the  saae  nuaber  of  degrees-of-freedom 
as  for  the  initial  analysis. 

The  points  for  each  contour  value  selected  were 
obtained  by  linearly  interpolating  between  the  nodal  values 
of  each  parameter  obtained  froa  the  initial  analysis.  The 
contours  were  constructed  by  saoothly  connecting  the  points 
by  hand.  The  element  layout  along  the  contours  posed  the 
aost  formidable  problem  because  the  coarse-to-fine  tran¬ 
sition  often  resulted  in  severe  eleaent  distortion,  and  it 
sometimes  became  necessary  to  degenerate  quadrilateral 
elements  into  triangles  when  the  transition  was  acute.  It 
was  decided  that  the  optimal  element  shapes  should  be 
preserved  along  the  contours  in  regions  of  highest  strass. 

The  contours  obtained  and  the  corresponding  grid  are 
presented  for  each  of  the  following  solution  resultants: 

•  torsion  function  (Fig.  5-6) 

•  shear  stress  (Fig.  6.7) 

•  strain  energy  density,  SED  (Fig.  6.8) 

The  resulting  grid  for  each  of  the  response  function 
contours  produces  smaller  elements  in  the  region  of  greatest 
stress  near  the  bottom  of  the  kayway  and  around  the 
periphery  of  the  shaft  where  the  stress  is  moderately  high. 
Consequently,  the  elements  near  the  center  where  the  stress 
is  zero  are  larger.  These,  of  course,  are  the  desired 
effects  for  an  optimization  criterion.  A  somewhat  unusual 
behavior  is  observed  at  the  point  of  intersection  of  the 
keyway  and  the  shaft  boundary  where  the  stress  is  also  zero. 
Apparently,  the  shear  stress  gradient  is  larger  than  the 
gradients  in  torsion  function  and  the  strain  energy  density. 


(b)  Corresponding  Grid 


Contouring  for  tbe  SED  Funct 


resulting  in  smaller  elements  being  produced  in  that  region 
by  the  shear  stress  criterion  (Pig.  6.7)  than  by  the  other 
criteria.  While  all  of  the  grids  possess  to  soae  extent  the 
desirable  features  cf  an  optimal  grid,  the  strain  energy 
density  function  produces  a  far  sore  drastic  refineaent 
towards  the  point  of  aaxiaua  stress,  while  the  others  repre- 
sent  aore  moderate  refinements.  In  fact,  the  SED  contours 
are  sc  dense  around  the  keyway  that  the  coarse-to-f ine 
element  transition  scheme  must  include  degenerate  quadrilat- 
erals  to  avoid  violating  the  contours.  Note  also  that  the 
coarse- tc-f ine  transition  for  the  torsion  function  response 
is  fairly  smooth  whereas  the  transition  for  the  strain 
energy  density  and  shear  stress  refinement  tends  to  produce 
distorted  and  elongated  elements.  This  is  aggravated  by  the 
fact  that,  unlike  the  torsion  function,  the  shear  stress  and 
strain  energy  density  are  not  monotonic  over  the  domain. 

The  solution  results  obtained  for  each  grid  are 
presented  in  the  upper  half  of  Table  IV.  At  first  glance, 
the  results  of  the  refinements  are  disappointing  in  compar¬ 
ison  to  the  parenthetic  values  obtained  using  a  uniform 
grid.  While  all  three  criteria  produce  improved  accuracy  for 
the  maximum  shear  stress,  the  errors  for  the  maximum  torsion 
function  value  grow  extremely  large.  Recalling  the  observa¬ 
tions  made  for  the  cne-diaensional  study,  it  would  follow 
that  such  behavior  is  probably  due  to  the  unusually  large 
elements  near  the  center  of  the  domain.  The  entries  in  the 
lower  half  of  Table  IV  reflect  the  drastic  improvement 
obtained  by  simply  introducing  a  few  additional  degrees-of- 
freedom  along  those  element  edges  which  grew  exceptionally 
long  during  the  optimization  process,  thus  increasing  their 
polynomial  order  from  one  to  two.  Wot  only  did  this  modifi¬ 
cation  reduce  the  large  errors  for  the  maximum  torsion 
function  estimation,  but  modest  improvements  were  3 Iso 
obtained  in  the  estimations  of  the  other  resultants  as  well. 


TABLE  IV 


‘Parenthetic  values  for  unifora  linear  grid  of  sane  nuaber  of 
degrees-of-f  reedca. 


Once  again,  the  selection  of  the  optimum  grid  would  depend 
predominantly  on  the  optiaization  goal  of  the  analysis,  bur 
one  would  likely  agree  that  the  strain  energy  density  varia¬ 
tion  along  with  some  modification  to  restrict  excessive 
eleaent  growth  provides  the  superior  refinement  criterion. 

in  additional  word  of  caution  is  in  order  for  the 
contouring  techniques  for  grid  optimization.  Because  the 
problea  aust  be  completely  redefined  froa  scratch  after  the 
initial  analysis,  the  preprocessing  cost  can  become  enor¬ 
mous,  especially  if  several  cycles  are  employed  to  obtain 
aore  precise  contours  as  some  authors  suggest.  Unless  there 
is  available  an  interactive  automatic  aesh  generator  based 
on  this  technique,  such  as  the  one  described  in  Reference 
[Hef.  11],  contouring  should  be  abandoned  in  favor  of  some 
aore  easily  iapleaented  grid  optimization  techniques 
employing  similar  refinement  criteria. 

2 .  Selective  Refinement 

The  simplest  way  to  avoid  the  problems  encountered 
in  the  contouring  techniques  is  to  perform  the  initial  anal¬ 
ysis  on  a  reasonably  coarse  grid  and  then  to  selectively 
refine  those  elements  over  which  the  strain  energy  density 
varies  the  most.  The  critical  concern  then  arises  as  to  how 
coarse  the  initial  grid  should  be.  If  the  preprocessor 
employs  the  necessary  constraint  conditions  to  permit  the 
"father-tc- four- sons"  eleaent  subdivision  scheme  directly, 
or  if  hierarchical  refinement  is  employed,  then  the  initial 
grid  should  be  just  coarse  enouoh  to  adequately  define  the 
problem  and  to  limit  the  number  of  refinement  cycles  neces¬ 
sary.  The  latter  becomes  even  less  of  a  concern  if  iterative 
solvers  are  employed.  If,  on  the  other  hand,  coarse-tc-fine 
transition  schemes  are  used  to  implement  the  h-versicn  or 
only  lew  polynomial  order  elements  are  available  in  the 
p-version,  then  the  initial  grid  must  be  sufficiently  fine 


so  as  not  to  restrict  severely  the  aaount  of  refinement 
which  can  be  performed  in  any  given  cycle.  Unfortunately, 
the  conditions  under  which  this  investigation  was  conducted 
were  those  of  the  latter. 

a.  The  h-Version 

Selective  refinement  by  the  h-version  was 
performed  on  both  linear  and  quadratic  element  grids.  For 
the  linear  case,  the  initial  analysis  was  performed  on  a 
uniform  grid  of  55  elements,  72  aodes,  and  50  degrees-of- 
freedom.  The  initial  quadratic  analysis  employed  an  eight 
element  uniform  grid  of  37  nodes  and  20  degrees-of-freedom. 
The  reason  for  such  a  great  disparity  in  the  number  of 
elements  for  the  initial  analyses  is  that  subdividing  a 
quadratic  element  introduces  many  more  degrees-of-freedom 
than  the  subdivision  of  a  linear  element.  These  numbers  were 
chosen  to  provide  approximately  the  same  number  of  degrees- 
of-freedcm  for  the  optimum  grid  of  the  final  cycle  for  each 
case.  The  initial  analysis  is  performed  and  those  elements 
over  which  the  strain  energy  density  variation  is  signifi¬ 
cantly  greater  become  candidates  for  refinement.  The  refine¬ 
ment  is  performed  by  subdividing  each  candidate  element  into 
four  new  ones  by  constructing  a  coarse-to-fine  transition 
zone  of  "buffer"  elements  around  the  refined  region. 
Successive  analyses  and  selective  refinements  are  repeated 
until  the  maximum  element  strain  energy  density  variation  is 
approximately  that  of  the  remainder  of  the  grid.  The  process 
is  further  improved  when  the  nodal  values  of  the  strain 
energy  density  are  used  to  indicate  the  general  direction  in 
which  the  refinement  is  to  proceed.  This  permits  multiple 
refinements  in  the  same  cycle,  thereby  reducing  the  number 
of  cycles  required  tc  arrive  at  the  optimal  grid.  For  this 
problem,  the  linear  grid  required  two  refinement  cycles 
while  the  quadratic  grid  required  three  cycles.  The 
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(b)  Cycle  1 
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!  selective  refinement  process  is  depicted  in  Figures  6.9  and 

6.10  for  the  linear  and  quadratic  element  grid, 
respectively. 

The  solution  results  for  each  selective  refine- 

I 

|  sent  cycle  are  presented  in  Table  7.  The  most  impressive 

observation  to  be  made  is  the  significant  improvement  in  the 
maximum  shear  stress  estimate  for  successive  cycles.  While 
there  is  also  modest  improvement  in  the  accuracy  of  the 
torque  estimate  for  successive  cycles,  when  compared  to  the 
estimate  obtained  from  the  uniform  grid  of  the  same  number 
of  degrees-of-freedom  and  polynomial  order,  the  refinement 
estimate  cf  torque  is  slightly  poorar.  This  is  because  addi¬ 
tional  degrees-of-freedom  are  being  introduced  in  only  a 
small  region  of  the  domain  but  the  torque,  and  energy,  are 
global  quantities.  The  author  has  no  satisfactory  explana¬ 
tion  why  the  estimate  for  the  maximum  torsion  function 
improves  for  successive  refinements  of  the  linear  grid  but 
not  for  the  quadratic  case.  However,  as  has  already  been 
mentioned,  this  particular  solution  parameter  appears  very 
sensitive  to  such  prcblen  variables  as  nodal  placements  and 
element  shapes;  hence,  its  behavior  is  difficult  to  predict, 
even  when  the  refinement  is  applied  to  regions  remote  from 
the  point  where  the  maximum  torsion  function  value  occurs, 

|  as  was  the  case  in  these  examples.  For  computational 

reasons,  it  is  desirable  to  restrict  the  number  of  refine¬ 
ment  cycles  to  a  necessary  minimum.  In  this  example,  the 
quadratic  grid  required  an  additional  cycle  over  the  linear 
grid  but  thi3  is  because  it  is  necessary  to  perform  the 
initial  quadratic  analysis  using  far  fewer  degrees-of- 
freedom.  Therefore,  the  6arly  cycles  of  the  quadratic  anal¬ 
ysis  actually  represent  comparatively  smaller  problems. 


74 


A."’;.'  ;V,‘.  **«.  V*. 


M  +J 
■J  a 
a  a 

«  m 
H  • 
0 


1 

1 

1 

0 

1 

1 

o> 

r» 

1 

1 

• 

• 

1 

®  1 

1 

©  1 

1 

i 

1 

1 

1 

VI 

1 

P* 

VO 

r» 

1 

1 

VO 

H  » 

0 

CN 

CN 

1 

VO 

1 

• 

• 

• 

1 

O 

1 

CN 

CN 

CN 

1 

• 

l 

1 

1 

1 

1 

o 

I 

1 

1 

| 

1 

t 

1 

r 

• 

CN 

CN 

i 

©  1 

0 

Ov 

• 

a  i 

• 

• 

i 

i 

lO 

m 

i 

VI 

1 

i 

i 

i 

w- 

i 

M  1 

i 

(0  1 

r* 

© 

o 

i 

t- 

,0  1 

VO 

cn 

© 

i 

CO 

K  ' 

• 

• 

• 

i 

• 

i 

p* 

CN 

o 

i 

ov 

i 

i 

t 

1 

• 

i 

i 

1 

i 

i 

i 

i 

m 

© 

i 

i  i 

vO 

m 

i 

i 

o 

o 

i 

i 

• 

• 

i 

M  1 

o 

o 

i 

0  1 

-w 

i 

0  1 

i 

0 

CO 

CN 

CN 

cn 

cn 

i 

i 

29 

i 

o 

O 

o 

i 

• 

i 

• 

• 

> 

i 

© 

i 

o 

o 

o 

i 

i 

i 

1 

1 

i 

i 

•»<  i 
OO  I 

o 

VO 

m 

i 

i 

i 

© 

as  ca  i 

m 

VO 

(>■ 

i 

i 

i 

•0 

CN 

W  1 

i 

•H 

■TO  4 

CN 

CN 

0 

i 

M 

1- 

OO  t 

TO 

r» 

O' 

o 

i 

© 

cn 

XX  | 

tH 

i 

1-4 

i 

** 

u 

i 

a 

0  1 

i 

® 

•  0  1 

4J 

m 

CN 

i 

0 

© 

Or-j  • 

a 

m 

r« 

00 

i 

® 

ZU  1 

0 

» 

i—4 

0 

i 

W 

+>  1 

® 

i 

0  1 

9  I 

• 

i 

£ 

0  1 

•H 

i 

+> 

i—4 

S  ’ 

U 

0 

i 

0 

0 

e®  i 

0 

•H 

i 

u 

•H 

•rtH  1 

V 

•P 

CN 

• 

TO 

<MO  1 

a 

•H 

i 

0 

•H 

®  >1| 

•H 

a 

i 

3 

3 

B3U  1 

►4 

M 

i 

Oi 

4-4 

00  vO  CN 

CN  ^  *■ 

O  O  O 
•  •  • 
©  ©  © 

I  I  I 


o  o  o 


cn  m  0 

mom 
•  •  • 

p"  vo  m 

*  i  i 


*-  o 
*  i 


0  <M  «*■ 

O  O  O 

•  •  • 

o  o  o 


«-  cn  in 


‘Parenthetic  values  for  unifora  grid  of  same  polynomial  order  and  same 
number  of  degcees-of -freedom. 


b.  The  p-Fersion 


Before  continuing  to  the  next  optimization  tech¬ 
nique,  it  is  worthwhile  to  take  a  quick  look  at  selective 
refineaent  employing  the  p-version.  Because  the  finite 
eleeent  code  used  in  this  investigation  only  provided  for 
element  orders  of  one  and  two,  the  advantages  of  the  method 
cannot  be  fully  realized,  but  the  effects  of  a  single  cycle 
can  be  examined. 

Beginning  with  three  uniform  grids  with 
differing  numbers  of  linear  elements,  the  initial  analyses 
were  performed.  In  each  case,  the  elements  over  which  the 
strain  energy  density  varied  the  most  were  transformed  from 
4-noded  Lagrangian  elements  into  8-noded  serendipity 
elements  by  the  addition  of  midsiie  nodes.  The  element 
grids  are  shown  in  Figure  6.11  and  the  asterisks  denote 
those  elements  for  which  the  polynomial  order  was  increased. 
In  this  example,  the  number  of  elements  to  be  refined  in 
each  case  was  chosen  so  as  to  achieve  approximately  the  same 
number  of  degrees-of -freedom  after  a  single  cycle. 

The  solution  results  are  shown  in  Table  ?i. 
Significant  improvements  in  the  estimate  of  the  maximum 
shear  stress  were  achieved  for  each  case.  An  improvement  in 
the  estimated  torque  was  also  realized  for  all  three  cases, 
but  was  more  noticeable  when  the  number  cf  refined  elements 
was  larger.  This  is  because  quadratic  elements  are  far 
superior  to  linear  elements  in  the  modeling  of  integral 
quantities,  as  observed  in  Figure  6.h.  Somewhat  disturbing 
is  the  increased  error  in  the  estimate  of  the  maximum 
torsion  function  value  observed  in  two  of  the  three  refine¬ 
ments  even  though  the  elements  in  the  vicinity  of  its  occur¬ 
rence  were  not  affected.  Again,  this  is  likely  attributable 
to  the  unusual  behavior  patterns  of  this  quantity  already 
discussed. 


(a)  76  Degraes-of- Fre  edoa  -  20  Refined  Elements 


(d)  76  Degrees- of -Freedom  -  9  Refined  Elements 


(c)  75  Oegrtes-of-Freedoa  -  4  Refined  Elements 


Figure  6.11  Selective  Refinement  Employing  p-Tersio 


•gas:  -aBMaftnagg.. 


‘Parenthetic  values  for  unifora  linear  grid  of  saae  number  of  degrees-of- 


Ia  order  to  perform  additional  cycles  of  the 
p-version,  it  would  be  necessary  to  alter  the  refinement 
criterion  slightly.  Because  the  element  sizes  do  not  change 
for  successive  cycles,  the  need  for  refinement  would  neces¬ 
sarily  be  based  on  strain  energy  density  variation  between 
nodes  rather  than  over  the  elements. 

Selective  refinement  employing  the  p-version  is 
most  efficiently  implemented  hierarchically,  in  which  case 
it  acquires  some  attractive  computational  advantages.  It  is 
unfortunate  that  time  did  not  permit  further  investigation 
here,  but  the  ne9d  for  future  research  is  evident. 

3 .  Subdsiajn  isolate 

The  refinement  criterion  and  initial  procedures  in 
employing  subdomain  isolation  are  identical  to  those  used  in 
selective  refinement.  After  the  candidate  elements  for 
refinement  are  identified,  they  are  completely  isolated  from 
the  remainder  of  the  domain  and  solved  as  smaller  subdomain 
problems.  The  advantages  of  the  technique  are  twofold.  By 
isolating  the  elements  to  be  refined  the  solution  is  not 
repeated  in  each  cycle  for  those  elements  for  which  the 
initial  analysis  is  assumed  adequate.  Furthermore,  by  elimi¬ 
nating  most  of  the  degrees-of-f reedom  over  the  entire 
domain,  the  subsequent  refinement  of  the  isolated  region  can 
be  much  greater  than  would  otherwise  be  practical. 

As  before,  the  technique  was  applied  to  both  linear 
and  quadratic  uniform  element  grids.  Those  elements  of  the 
initial  analysis  over  which  the  strain  energy  density  varia¬ 
tion  was  exceptionally  large  were  isolated  to  comprise  the 
subdomain  in  each  case.  There  were  three  such  elements  of 
the  initial  linear  grid  and  two  for  the  quadratic  grid.  Each 
subdomain  was  uniformly  refined  to  achieve  approximately  the 
same  number  of  degrees-of- freedom  as  the  initial  analyses. 
The  process  is  depicted  in  Figure  5.12  for  the  linear  grid 
and  Figure  6.13  for  the  quadratic  case. 


(b)  aefineasnt  of  Isolated  Subdomain 


Figure  6.12  Subdoaain  isolation  -  Linear  Eleaents 


(b)  Refinement  of  Isolated  Subdomain 


Pigure 


.13  Sabdoaain  Isolation  -  Quadratic 

8  1 


Eleaents. 


In  performing  subdomain  isolation  the  governing 
equations  remain  the  same,  while  only  the  domain  and  the 
boundary  conditions  are  altered.  When  the  subdomain 

boundary  has  nodes  common  to  the  initial  grid,  then  the 
boundary  values  for  those  nodes  are  simply  the  solution 
values  obtained  from  the  initial  analysis.  The  boundary 
values  arising  from  the  introduction  of  new  boundary  nodes 
during  the  refinement  process  must  be  generated  by  interpo¬ 
lation  of  the  soluticn  results  of  the  initial  analysis.  One 
of  the  options  for  an  interpolation  scheme  is  simply  to  use 
the  shape  functions  of  the  unrefined  elements.  This  may  not 
be  desirable  in  the  case  of  linear  elements,  so  a  higher 
order  interpolation  may  be  employed.  In  this  example  a  third 
order  Lagrangian  polynomial  was  convenient  for  the  linear 
case  since  there  are  four  nodes  from  the  initial  analysis 
along  the  right-hand  boundary  of  the  subdomain.  Since  there 
are  only  two  such  nodes  on  the  upper  boundary,  it  is  neces¬ 
sary  to  "borrow"  some  adjacent  nodal  values  from  the 
discarded  portion  of  the  domain  in  the  generation  of  new 
boundary  values  using  higher  order  interpolation. 

The  solution  results  presented  in  Table  VII  are 

quite  remarkable.  In  a  single  cycle,  the  solution  accuracy 
for  the  maximum  shear  stress  has  increased  by  a  full  order 
of  magnitude.  No  other  optimization  technique  examined  in 
this  investigation  produced  such  improvements.  Note  that 
the  higher  order  polynomial  interpolation  for  the  boundary 
values  did  improve  the  solution  results  for  the  linear  case. 
One  of  the  disadvantages  of  this  technique  is  that  the 

refinement  can  produce  no  improvement  in  the  estimation  of 
local  quantities  outside  the  subdomain.  As  in  this  example, 
the  estimation  of  the  maximum  torsion  function  value 

obtained  from  the  initial  analysis  must  be  accepted  as  the 

optimum  since  it  occurs  outside  the  subdomain.  Furthermore, 
since  its  value  is  predominantly  affected  by  refinements  in 
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regions  of  high  gradients,  it  is  doubtful  that  isolating  a 
nev  subdonain  using  the  elements  adjacent  to  its  point  of 
occurrence  would  noticeably  improve  its  accuracy.  However, 
since  the  torque  is  a  globally  computed  quantity,  refinement 
will  improve  the  accuracy  of  its  contribution  from  the 
subdomain  resulting  in  improvement  of  its  accuracy  overall 
as  observed  in  Table  VII.  It  is  this  strictly  local  nature 
of  the  subdomain  isolation  technique  which  restricts  its 
applicability.  But  if  the  optimization  goals  are  well 
defined  and  it  is  understood  under  which  conditions  and  for 
what  parameters  it  is  effective,  it  can  be  an  extremely 
powerful  grid  optimization  technique. 

4.  JL§sh  isadjaa 

While  mesh  grading  is  nearly  always  performed  on  an 
"a-priori"  basis,  it  may  also  be  employed  adaptively  to 
provide  a  simple  grid  optimization  technique.  After  an 
initial  solution  has  been  obtained,  the  analysis  may  be 
repeated  using  various  combinations  of  grading  ratios  in 
order  to  achieve  a  more  uniform  distribution  of  the  element 
strain  energy  density  variations.  Here  the  grading  ratio 
refers  to  the  constant  ratio  of  adjacent  element  lengths 
along  a  boundary  of  the  domain  to  which  grading  is  applied. 
There  are  several  drawbacks  to  the  technique,  the  first 
being  that  a  good  combination  of  grading  ratios  may  be 
difficult  to  obtain  in  a  reasonable  number  of  cycles.  The 
other  difficulty  is  that  if  smaller  elements  are  desired  in 
more  than  one  region  the  domain  must  be  refined  and 
constructed  by  subregions. 

Unfortunately,  the  domain  of  this  problem  is  not 
well  suited  for  mesh  grading  since  it  possesses  no  favorable 
geometric  symmetry.  Hence,  the  resulting  element  elongation 
and  distortion  would  become  severe  for  larger  grading 
ratios.  For  simplicity,  the  nodal  placements  will  be  biased 
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only  along  tha  tvo  domain  boundaries  adjacent  to  the  point 
of  maximum  stress  and  the  same  grading  ratio  will  be  applied 
for  each.  This  will  result  in  snail  eleaents  near  the  bottom 
of  the  keyvay  and  large  elements  along  the  periphery  of  the 
shaft.  Hhile  this  is  not  the  most  desirable  grid  topology, 
it  will  produce  a  aore  uniform  distribution  of  the  element 
SBD  variations. 

The  technique  was  applied  to  both  linear  and  quad¬ 
ratic  element  grids  starting  with  a  uniform  mesh  and  succes¬ 
sively  increasing  the  grading  ratio  until  the  elements  along 
the  shaft  periphery  exhibited  SED  variations  as  large  as 
those  for  the  eleaents  along  the  keyvay.  In  both  cases  this 
condition  occurred  beyond  the  point  where  excessive  element 
elongation  would  be  expected  to  produce  numerical  inaccura¬ 
cies.  Graded  meshes  for  selected  grading  ratios  r,  are  shown 
in  Figure  6.1U  for  linear  elements  and  Pigure  6.15  for  quad¬ 
ratic  elements.  The  solution  results  are  presented  in 
Table  VIII.  As  can  be  observed,  the  maximum  shear  stress 
estimate  improves  for  each  successive  increase  in  the 
grading  ratio.  However,  the  cost  of  such  improvement  is  the 
accompanying  degradation  in  the  estimate  of  the  maximum 
torsion  function  value.  This  is  to  be  expected  since  the  tvo 
maxima  occur  at  different  locations  in  the  domain  and  there¬ 
fore  decreasing  the  size  of  the  eleaents  in  the  vicinity  of 
one  will  necessarily  increase  the  element  size  near  the 
other.  Hote  that  this  degradation  is  not  nearly  so  severe 
for  the  quadratic  eleaents,  and  the  accuracy  actually 
improves  for  a  low  value  of  the  grading  ratio.  This  is 
because  the  higher  order  interpolation  can  better  accomodate 
larger  elements.  For  both  linear  and  quadratic  element  grids 
the  optimal  accuracy  in  torque  estimation  occurs  for  the 
moderately  graded  meshes. 


Figure  6.14  Mesh  Grading  -  Linear  Elements 


Figure  6.15  Hesfa  Grading  -  Quadratic  Elements 


TABLE  VIII 

flesh  Grading 

Solution  Besults 

Pare 

entage  Relative  Er 

ror 

Grading 

Hatio 

^aa  x 

Toax  /G9 

I  /  G0 

Linear  Element  Grids 

(  78  elements,  98 
72  degrees-of-f 

nodes, 
reedom  ) 

1.0 

0.060 

-6.  06 

-1.77 

1.1 

0.161 

-2.63 

-1.56 

1.2 

0.389 

-1.03 

-1.77 

1.3 

0.679 

-0.  484 

-2.  20 

Quadratic 

Bleaent  Grids  (  28  elenents,  107  nodes, 

76  degrees-of-freedom  ) 

1.0 

-0.0093 

-5.26 

-0.0116 

1.1 

0.0063 

-1.85 

-0.0064 

1.2 

0.0162 

-0.606 

-0. 0091 

1.3 

-0.0246 

-0.413 

-0.0179 

It  is  likely  that  some  of  the  error  is  attributable 
to  the  numerical  inaccuracies  due  to  eleaent  distortion, 
ihen  applying  a  grading  technique  the  analyst  should  seek  an 
equitable  balance  between  the  refineaent  criterion  and  the 
grid  topclogy. 
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Finally,  since  the  grading  ratio  is  usually  applied 
to  the  nodal  separation  rather  than  the  element  edge 
lengths,  it  is  advisable  to  reposition  the  aidside  nodes  of 
guadratic  elements  so  that  they  lie  near  the  center  of  the 
element  edges.  This  will  generally  improve  the  accuracy  of 
all  the  solution  parameters,  especially  if  the  grading  is 
somewhat  extreme. 

2.  GUIDELINES  FOR  T  BO-DIHE  US  I  ORAL  GRID  OPTIHIZATIOH 

In  order  to  provide  some  guidelines  for  obtaining 
optimal  finite  element  solutions  for  two-dimensional  prob¬ 
lems  it  is  helpful  tc  compare  the  solution  results  obtained 
for  this  problem  employing  the  optimization  techniques 
available.  Such  a  comparison  is  presented  in  Table  IX.  The 
upper  portion  is  for  those  techniques  for  which  the  initial 
analysis  was  performed  using  linear  elements  and  the  lower 
portion  using  quadratic  elements.  Rote  that  all  of  the  grids 
employ  approximately  the  same  number  of  degrees-of- freedom, 
which  was  the  chosen  measure  of  analysis  cost  in  this 
investigation. 

In  making  this  comparison  it  is  important  to  understand 
just  how  significant  a  change  in  error  actually  is.  If  the 
convergence  rate  of  a  solution  parameter  is  very  slow,  even 
a  small  reduction  in  the  error  may  require  many  more 
degr 6 es-cf- freedom.  For  this  reason,  the  numbers  in  paren¬ 
theses  have  been  included  by  each  error  to  provide  a  rough 
approximation  of  the  number  of  degrees-of-freedom  required 
to  obtain  similar  accuracy  using  a  uniform  grid  of  the  same 
number  of  degrees-of-freedom  and  elements  of  the  same  poly¬ 
nomial  order  as  the  initial  analysis.  In  some  cases,  the 
analyses  were  not  actually  performed  but  instead  the  numbers 
in  parentheses  were  obtained  by  extrapolation  of  the  error 
versus  degrees-of-freedom  curve  (Fig.  6.4),  and  ignoring 
round-off  and  ill-conditioning  effects. 
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‘Parenthetic  values  indicate  the  approximate  number  of  degrees-of- freedom 
red  to  achieve  similar  accuracy  using  a  uniform  grid  of  the  same  order  as 
nitial  analysis. 


The  first  observation  to  be  made  from  Table  IX  is  that 
while  all  of  the  optimization  techniques  produced  signifi¬ 
cant  improvements  in  the  accuracy  of  the  maximum  derivative 
quantity  Tmax,  the  same  cannot  be  said  for  the  maximum 
solution  quantity  l/^max  and  the  integral  quantity  T.  One 
might  even  conclude  that  grid  optimization  is  not  cost 
effective  in  the  computation  of  these  values  since  the 
uniform  grid  provides  estimates  which  are  nearly  as  accu¬ 
rate,  and  in  some  cases  better,  than  the  optimum  grid.  This 
conclusion  would  be  correct  if  the  solution  maximum  and  its 
integral  were  the  only  resultants  of  interest  in  performing 
the  analysis.  Since  the  purpose  of  this  study  was  to  find  an 
optimum  grid  which  produced  acceptable  errors  for  all  the 
resultants,  the  uniform  grid  is  clearly  inadequate. 
Moreover,  since  in  the  majority  of  engineering  problems  it 
is  the  derivative  cf  the  solution  variable  which  is  of 
primary  interest,  it  deserves  special  consideration  in 
making  this  comparison. 

Furthermore,  one  might  conclude  that  the  reason  the 
error  in  the  maximum  solution  variable  is  larger  for  the 
optimal  grid  is  because  the  strain  energy  density  variation 
criterion  always  concentrates  the  degrees-of -freedom  in  the 
vicinity  of  the  maximum  derivative  value,  which  in  this  case 
does  not  coincide  with  that  of  the  maximum  solution  variable 
value.  However,  such  a  conclusion  is  incorrect  and  might 
erroneously  lead  one  to  attempt  refinement  where  the  maximum 
solution  variable  occurs  in  an  effort  to  improve  its  esti¬ 
mate.  Other  investigations  have  revealed  that  in  nearly  all 
cases  the  maximum  accuracies  for  all  three  solution  resul¬ 
tants  are  obtained  by  refinement  in  the  regions  of  highest 
strain  energy  density  variation  [lef.  11].  The  more  likely 
source  of  increasing  error  for  the  optimal  grids  is  the 
element  distortion  which  was  encountered  in  all  but  two 
techniques,  selective  p-version  refinement  and  subdomain 
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isolation.  Such  distortion  can  be  avoided  but  would  require 
aore  sophisticated  refinesent  techniques  than  were  available 
for  this  investigaticn. 

A  reasonable  choice  for  the  optiaun  grid  in  Table  IX 
would  be  one  for  which  all  three  values  in  parentheses  are 
as  large  as  reasonably  possible  taking  into  consideration 
the  nuaber  of  cycles  required  to  provide  such  accuracy. 
Based  on  such  a  criterion,  the  author  is  partial  to  subdo- 
aain  isolation  for  the  solution  of  two-dimensional  probleas 
using  linear  elements,  and  selective  refinement  for  finite 
eleaent  solutions  using  quadratic  eleaents.  clearly,  before 
a  concrete  recoaaendation  could  be  aade  for  a  wide  range  of 
applications,  aany  aore  problems  would  have  to  be  studied, 
but  these  two  techniques  were  fairly  simple  to  implement  for 
a  standard  finite  element  code  and  the  accelerated  conver¬ 
gence  of  the  solution  resultants  of  interest  was  impressive. 
Conceivably,  even  greater  solution  accuracies  might  be 
obtained  by  using  two  or  aore  of  these  techniques  in 
coabination. 

Here  again  the  crucial  eleaent  in  selecting  the  proper 
optimization  strategy  is  the  precise  definition  of  the 
purpose  for  which  the  finite  eleaent  analysis  is  to  be 
performed.  The  results  of  Table  IX  tend  to  support  the 
following  recoaaendations  and  conclusions: 

(1)  Regardless  of  the  optiaization  strategy  chosen, 
higher  order  eleaents  are  indispensable  if  high 
accuracies  for  integral  solution  quantities  are 
desired. 

(2)  If  the  aaziaua  derivative  of  the  solution  variable 
is  of  greatest  concern,  the  strictly  local  refine- 
eents  employing  subdoaain  isolation  techniques  can 
provide  exceptional  accuracy  for  a  minimum  number  of 
refineaent  cycles. 


(3)  If  the  maximum  solution  variable  value  occurs  at  a 
point  in  the  domain  removed  from  the  vicinity  of  the 
maximum  derivative  value,  then  its  best  estimate 
will  likely  be  obtained  using  a  reasonably  fine 
uniform  grid  and  selectively  subdividing  elements  in 
the  regions  cf  large  strain  energy  density 
variations. 
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The  purpose  of  this  paper  has  been  to  present  an  over¬ 
view  of  some  readily  eaployed  finite  element  grid  optimiza¬ 
tion  methods  and  to  demonstrate  their  effectiveness  in  the 
application  to  some  siaple  probleas.  This  work  is  by  no 
aeans  all  inclusive  and  the  subject  is  still  in  its  infancy. 
While  there  are  many  competing  approaches  to  the  problem, 
there  is  auch  aore  research  to  be  done  before  any  one 
becomes  widely  accepted  as  a  standard  analysis  tool.  Because 
of  the  limited  tiae  and  resources  available,  some  of  the 
aore  sophisticated  refinement  criteria  and  techniques  which 
have  been  developed  have  not  been  examined  in  any  detail. 
Instead,  the  approach  has  been  to  examine  those  techniques 
which  can  be  easily  incorporated  in  a  basic  finite  element 
coda.  However,  it  is  likely  that  some  recently  developed  and 
rather  elaborate  self-adaptive  grid  optimization  codes  will 
soon  be  available. 

Also,  this  paper  has  not  rained  the  important  classes 
of  probleas  in  dynaaic  and  nonlinear  analysis.  There  is 
considerable  ongoing  research  in  the  extension  of  these 
techniques  to  such  problems,  but  the  increased  complexity  is 
evident.  For  example,  in  vibration  analysis  there  is  an 
optiaua  grid  for  each  unique  eigenvalue,  but  it  is  for  these 
types  of  prcbleas  that  grid  optimization  is  most  promising. 

At  the  beginning  of  this  paper  it  was  stated  that  the 
goa.1  of  grid  optimization  was  to  obtain  maximum  solution 
accuracy  for  a  given  analysis  cost.  Throughout  this  paper  it 
has  been  shown  that,  prior  to  successfully  embarking  upon 
such  a  strategy,  the  underlying  purpose  of  the  analysis  must 
be  explicitly  defined.  Hopefully,  it  has  been  demonstrated 
that  grid  optimization  is  by  no  means  an  unrealistic  goal 


and  is  far  more  attractive  than  the  non-adaptive  practices 
widely  used  today. 

The  following  are  reco emendations  for  future  research 
topics: 

(1)  Investigation  of  aore  sophisticated  refinement 
criteria  based  on  eleaent  residuals  and  reliable 
error  estimates. 

(2)  Investigation  of  grid  optimization  techniques 
employing  adaptive  application  of  the  p-version. 

<3 )  Implementation  of  a  finite  element  preprocessor  for 
performing  hierarchical  grid  refinement. 

(4)  Implementation  of  a  self-adaptive  finite  element 
code. 


(5)  Application  of  grid  optimization  techniques  to  prob¬ 
lems  of  dynamic  and  nonlinear  analysis. 
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PORHULATIOH  OF  THE  ELASTIC  CABLE  PROBLEH 

A.  PROBLEM  ST  AT  EBERT 

Consider  a  perfectly  elastic  cable  initially  stretched 
between  two  fixed  points  a  distance  2L  apart  and  under 
tension  T.  If  the  cable  bears  a  distributed  load  per  unit 


Figure  A. 1  Tension  Cable  Under  Distributed  Loading. 

length  f  (x)  a3  shown  in  Figure  A.1,  the  governing  differen 
tial  equation  for  the  downward  defection  v(x)  is: 

2 

T  ♦  f(x)  »  0  (Eqn.  A.1) 


subject  to  the  essential  boundary  condition: 


?  (x  »  ±L)  «  0 


(Egn.  A  .2) 


If  the  distributed  load  is  a  supportive  load  provided  by 
a  iinkler  foundation  cf  nod  ulus  k  such  that 

f(x)  »  -  k  v(x) 

and  if  a  concentrated  load  P  is  applied  at  the  oidspan. 
Equation  A.  1  becoaes: 


(Eqn.  A. 3) 


subject  to  the  natural  boundary  condition: 


_  P/2 
T 


(Eqn.  A. 4) 


B.  PBOBLEH  SO  LOTI  OH 

The  analytical  solution  of  the  two-point  boundary  value 
prcblea  is: 


v  (x) 


C  tanh  X  L  cosh  X  x  -  sinh  X  x] 


where  X 


(0  £  x  <  L) 
(Eqn.  A. 5) 


The  finite  eleaent  solution  is  obtained  by  the  Galerkin 
foraulation  using  the  consistent  rather  than  the  lumped 
approxiaation  for  the  distributed  loading. 
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FOBBOLAIXOB  07  TIE  TAPBBBD  BAB  PBOBLZfl 


A.  PBOBLEH  STAIEHEBT 

Consider  a  tapered  bar  of  length  L  and  constant  aodulns 
of  elasticity  E  fixed  at  one  end.  The  cross-sectional  area 
A(x)  varies  linearly  froa  Aq  at  the  fixed  end  to  &t  at  the 
tip.  Let  the  bar  be  loaded  axially  by  a  concentrated  tip 
load  P#  and  a  distributed  load  for  which  the  intensity  g  (x) 
varies  linearly  froa  g  at  the  fixed  end  to  g  at  the  tip  as 


Flgere  B.l  Tapered  Bar  vith  Applied  Loads. 


shewn  in  figure  B.l.  The  governing  differential  agnation 
for  the  axial  displaceaent  u(x)  is: 


BjfefhStel  ♦  q  ■  o 


(0  S  X  5  L) 


(Egn.  B.l) 


subject  to  the  essential  boundary  condition: 


u  (x  *  0)  »  0 


(Egn.  B.2) 


and  the  natural  boundary  condition: 


(Egn.  B.3) 


B.  PBOBLBB  SOLOTIOI 
Lot 


«  *  1  ’  *t 


and 


i  -  qt/q 


o* 


u(x) 


For  P(x)  ■  P  ♦ 

X  £00 

o  Vx> 


if; 


£ 


dx 


g(x)  dx  tha  solution  is: 


aA  E 
o 


2§U-  Jx>  -  L(1  o 


1*  P  ax.  .  ,,  0  .  gx" 

~  ~  ~f)  +  (1-  ~ )*  -  ~ 

L 


2a'"  4L 
(0  <  x  <  L) 1 


(Egn.  B.4) 


Tbs  finits  slsssnt  solution  is  obtained  by  the  Galerkin 
formulation  using  the  consistent  rather  than  the  lusped 
approxiaation  for  the  distributed  loading. 


FOBBDLATIOB  OF  IBB  TOBSZOB  PBOBLBB 

i.  PBOBLBB  SI  AT  EH  BIT 

Consider  a  solid  circal&r  shaft  of  radius  "a"  and  shear 
aodulus  G,  with  a  circular  groove,  or  kayvay,  of  radius  b 
along  a  generator  of  the  shaft  with  the  cross-section  shown 


y 


Figure  C. 1  cross-section  of  Shaft  with  Keyway. 

in  Figure  c.1.  An  applied  torgue  I  will  produce  an  angle  of 
twist  per  unit  length  9.  The  aguilibriuw  condition  is 
satisfied  if  a  torsional  stress  function  $  exists  such  that 
the  shear  stress  components  are: 

t  *  SB  1^-  and  T  »  -G6  r^- 


The  governing  differential  equation  for  the  torsional  stress 
function  i>  ( x,j )  is: 

lit  +  +  2  *  0  (Eqn  •  C.1) 

3x2  3y2 

subject  to  the  Oirichlet  condition*  $  *  0  on  the  boundary. 

B.  PIOBIBH  SOLUTIOB 

The  solution  of  Equation  C.1  [Ref.  22:  pp.  141-143]  is: 

2 

4t(x,y)  ■  a(x-b2 — * —  )  -  h(x2+y2)  +  “  (Eqn.  C.2) 

x2+y2 

The  uazisue  shear  stress  ocurs  at  point  k  (Fig.  C.1)  and  is: 


Tuaz  *  G8  (2a  -  b) 


(Eqn.  C.3) 


The  applied  torque  conputed  from  the  area  integral: 


L 


T  *  2G0  /  ^  d k 
A 


g«  r 

4  L 


(4a4-8a2b2-2b4)a  +  (2a3b+7ab3  )4-tn 


■] 


( Eqn  C.  4) 


where  a  «  arcos  (b/2a) . 

The  strain  energy  per  unit  length  of  shaft  is: 


"  *  kfT* 


<JA  »  T9 


(Eqn.  C.5) 
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The  variational  foraulation  of  the  finite  eleoent 
approxiaation  is  presented  in  detail  in  Chapter  6  of 

Bafaraace  25. 
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